
#include "../include/OptionUtil.h"

double OptionUtil::ND(double x)
{
    return 1.0 / sqrt(2.0 * 3.141592653590) * exp(-0.5 * x * x);
}

double OptionUtil::CND(double x)
{
    double y, Exponential, SumA, SumB;
    double result;
    y = fabs(x);
    if (y > 37.0)
    {
        return 0.0;
    }
    else
    {
        Exponential = exp(-0.5 * y * y);
        if (y < 7.07106781186547)
        {
            SumA = 3.52624965998911E-02 * y + 0.700383064443688;
            SumA = SumA * y + 6.37396220353165;
            SumA = SumA * y + 33.912866078383;
            SumA = SumA * y + 112.079291497871;
            SumA = SumA * y + 221.213596169931;
            SumA = SumA * y + 220.206867912376;
            SumB = 8.83883476483184E-02 * y + 1.75566716318264;
            SumB = SumB * y + 16.064177579207;
            SumB = SumB * y + 86.7807322029461;
            SumB = SumB * y + 296.564248779674;
            SumB = SumB * y + 637.333633378831;
            SumB = SumB * y + 793.826512519948;
            SumB = SumB * y + 440.413735824752;
            result = Exponential * SumA / SumB;
        }
        else
        {
            SumA = y + 0.65;
            SumA = y + 4.0 / SumA;
            SumA = y + 3.0 / SumA;
            SumA = y + 2.0 / SumA;
            SumA = y + 1.0 / SumA;
            result = Exponential / (SumA * 2.506628274631);
        }

        if (x > 0) result = 1.0 - result;
        return result;

    }

}

double OptionUtil::CND2(double x)
{
    if (x == 0.0)
    {
        return 0.5;
    }
    else
    {
        double L, k;
        double result;
        double a1 = 0.31938153;
        double a2 = -0.356563782;
        double a3 = 1.781477937;
        double a4 = -1.821255978;
        double a5 = 1.330274429;
        L = fabs(x);
        k = 1.0 / (1.0 + 0.2316419 * L);
        result = 1.0 - 1.0 / sqrt(2.0 * 3.141592653590) * exp(-0.5 * L * L) * (a1 * k + a2 * k * k + a3 * k * k * k + a4 * k * k * k * k + a5 * k * k * k * k * k);
        if (x < 0.0) result = 1.0 - result;
        return result;

    }
}

double OptionUtil::BND(double X, double y, double rho)
{
    return 1.0 / (2.0 * 3.141592653590 * sqrt(1.0 - rho * rho)) * exp(-1.0 / (2.0 * (1.0 - rho * rho)) * (X * X + y * y - 2.0 * X * y * rho));
}

double OptionUtil::CBND(double X, double y, double rho)
{
    //'
    //'     A function for computing bivariate normal probabilities.
    //'
    //'       Alan Genz
    //'       Department of Mathematics
    //'       Washington State University
    //'       Pullman, WA 99164-3113
    //'       Email : alangenz@wsu.edu
    //'
    //'    This function is based on the method described by
    //'        Drezner, Z and G.O. Wesolowsky, (1990),
    //'        On the computation of the bivariate normal integral,
    //'        Journal of Statist. Comput. Simul. 35, pp. 101-107,
    //'    with major modifications for double precision, and for |R| close to 1.
    //'   This code was originally transelated into VBA by Graeme West
    int i, ISs, LG, NG;
    //Dim XX(10, 3) As Double, W(10, 3) As Double
    double XX[11][4];
    double W[11][4];

    double h, k, hk, hs, BVN, Ass, asr, sn;
    double A, b, bs, c, d;
    double xs, rs;


    W[1][1] = 0.17132449237917;
    XX[1][1] = -0.932469514203152;
    W[2][1] = 0.360761573048138;
    XX[2][1] = -0.661209386466265;
    W[3][1] = 0.46791393457269;
    XX[3][1] = -0.238619186083197;

    W[1][2] = 4.71753363865118E-02;
    XX[1][2] = -0.981560634246719;
    W[2][2] = 0.106939325995318;
    XX[2][2] = -0.904117256370475;
    W[3][2] = 0.160078328543346;
    XX[3][2] = -0.769902674194305;
    W[4][2] = 0.203167426723066;
    XX[4][2] = -0.587317954286617;
    W[5][2] = 0.233492536538355;
    XX[5][2] = -0.36783149899818;
    W[6][2] = 0.249147045813403;
    XX[6][2] = -0.125233408511469;

    W[1][3] = 1.76140071391521E-02;
    XX[1][3] = -0.993128599185095;
    W[2][3] = 4.06014298003869E-02;
    XX[2][3] = -0.963971927277914;
    W[3][3] = 6.26720483341091E-02;
    XX[3][3] = -0.912234428251326;
    W[4][3] = 8.32767415767048E-02;
    XX[4][3] = -0.839116971822219;
    W[5][3] = 0.10193011981724;
    XX[5][3] = -0.746331906460151;
    W[6][3] = 0.118194531961518;
    XX[6][3] = -0.636053680726515;
    W[7][3] = 0.131688638449177;
    XX[7][3] = -0.510867001950827;
    W[8][3] = 0.142096109318382;
    XX[8][3] = -0.37370608871542;
    W[9][3] = 0.149172986472604;
    XX[9][3] = -0.227785851141645;
    W[10][3] = 0.152753387130726;
    XX[10][3] = -7.65265211334973E-02;

    if (fabs(rho) < 0.3)
    {
      NG = 1;
      LG = 3;
    }
    else if( fabs(rho) < 0.75)
    {
      NG = 2;
      LG = 6;
    }
    else
    {
      NG = 3;
      LG = 10;
    }


    h = -X;
    k = -y;
    hk = h * k;
    BVN = 0.0;

    if (fabs(rho) < 0.925)
    {
        if (fabs(rho) > 0.0)
        {
            hs = (h * h + k * k) / 2.0;
            asr = ArcSin(rho);
            for (i = 1; i <= LG; i++)
            {
                for (ISs = -1; ISs <= 1; i = i + 2)
                {
                    sn = sin(asr * (ISs * XX[i][NG] + 1) / 2.0);
                    BVN = BVN + W[i][NG] * exp((sn * hk - hs) / (1.0 - sn * sn));
                }
            }
            BVN = BVN * asr / (4.0 * 3.141592653590);
        }
        BVN = BVN + CND(-h) * CND(-k);
    }
    else
    {
        if (rho < 0.0)
        {
            k = -k;
            hk = -hk;
        }
        if (fabs(rho) < 1.0)
        {
            Ass = (1.0 - rho) * (1.0 + rho);
            A = sqrt(Ass);
            bs = (h - k) *(h - k);
            c = (4 - hk) / 8.0;
            d = (12 - hk) / 16.0;
            asr = -1.0*(bs / Ass + hk) / 2.0;
            if (asr > -100.0) BVN = A * exp(asr) * (1.0 - c * (bs - Ass) * (1.0 - d * bs / 5) / 3.0 + c * d * Ass * Ass / 5.0);
            if (-hk < 100)
            {
              b = sqrt(bs);
              BVN = BVN - exp(-hk / 2) * sqrt(2.0 * 3.141592653590) * CND(-b / A) * b * (1.0 - c * bs * (1.0 - d * bs / 5.0) / 3.0);
            }
            A = A / 2.0;
            for(i = 1; i<= LG; i++)
            {
                for( ISs = -1 ; ISs<=1; ISs = ISs+2)
                {
                    xs = (A * (ISs * XX[i][NG] + 1.0)) * (A * (ISs * XX[i][NG] + 1.0));
                    rs = sqrt(1.0 - xs);
                    asr = -(bs / xs + hk) / 2.0;
                    if( asr > -100)
                    {
                        BVN = BVN + A * W[i][NG] * exp(asr) * (exp(-hk * (1.0 - rs) / (2.0 * (1.0 + rs))) / rs - (1.0 + c * xs * (1.0 + d * xs)));
                    }
                }
            }
            BVN = -BVN / (2 * 3.141592653590);
        }

        if (rho > 0.0)
        {
            BVN = BVN + CND(-Max(h, k));
        }
        else
        {
            BVN = -BVN;
            if (k > h) BVN = BVN + CND(k) - CND(h);
        }
    }


    return BVN; //CBND = BVN
}

double OptionUtil::CBND2(double A, double b, double rho)
{
    double g, P, sum;

    double X[5] = { 0.018854042, 0.038088059, 0.0452707394, 0.038088059, 0.018854042};
    double y[5] =  { 0.04691008, 0.23076534, 0.5, 0.76923466, 0.95308992 };
    sum = 0.0;
    for (int i = 0; i < 5; i++)
    {
       P = y[i] * rho;
       g = 1 - P * P;
       sum = sum + X[i] * exp(0.5 * (2.0 * A * b * P - A * A - b * b) / g) / sqrt(g);

    }

    return rho * sum + CND(A) * CND(b);

}

double OptionUtil::CBND3(double A, double b, double rho)
{
    double X[5] =  { 0.24840615, 0.39233107, 0.21141819, 0.03324666, 0.00082485334};
    double y[5] =  { 0.10024215, 0.48281397, 1.0609498, 1.7797294, 2.6697604 };
    double rho1, rho2, delta;
    double a1, b1, sum;
    int i, j;
    a1 = A / sqrt(2.0 * (1.0 - rho * rho));
    b1 = b / sqrt(2.0 * (1.0 - rho * rho));
    double result=0.0;

    if (A <= 0.0 && b <= 0.0 && rho <= 0.0)
    {
        sum = 0.0;
        for( i = 0; i<5; i++)
        {
            for(j=0; j<5; j++)
            {
                sum = sum + X[i] * X[j] * exp(a1 * (2.0 * y[i] - a1)
                    + b1 * (2.0 * y[j] - b1) + 2.0 * rho * (y[i] - a1) * (y[j] - b1));
                result = sqrt(1.0 - rho * rho) / 3.141592653590 * sum;
            }
        }

    }
    else if(A <= 0.0 && b >= 0.0 && rho >= 0.0)
    {
        result = CND(A) - CBND3(A, -b, -rho);
    }
    else if (A >= 0.0 && b <= 0.0 && rho >= 0.0)
    {
        result = CND(b) - CBND3(-A, b, -rho);
    }
    else if (A >= 0.0 && b >= 0.0 && rho <= 0.0)
    {
        result = CND(A) + CND(b) - 1 + CBND3(-A, -b, rho);
    }
    else if (A * b * rho > 0.0)
    {
        rho1 = (rho * A - b) * Sign(A) / sqrt(A * A - 2 * rho * A * b + b * b);
        rho2 = (rho * b - A) * Sign(b) / sqrt(A * A - 2 * rho * A * b + b * b);
        delta = (1 - Sign(A) * Sign(b)) / 4.0;
        result = CBND3(A, 0.0, rho1) + CBND3(b, 0.0, rho2) - delta;
    }


    return result;
}

double OptionUtil::CBND4(double A, double b, double rho)
{
    //'modified/corrected from the second function in Drez & Wes paper pg. 105
            //'0/0 case resolved by l'H rule

    int i;

    double X[5] =  { 0.04691008, 0.23076534, 0.5, 0.76923466, 0.95308992 };
    double W[5] =  { 0.018854042, 0.038088059, 0.0452707394, 0.038088059, 0.018854042 };

    double h1, h2;
    double LH=0.0, h12, h3, h5, h6, h7, h8;
    double r1, r2, r3, rr;
    double AA, ab;

    h1 = A;
    h2 = b;
    h12 = (h1 * h1 + h2 * h2) / 2.0;

    double result = 0.0;

    if (fabs(rho) >= 0.7)
    {
        r2 = 1.0 - rho * rho;
        r3 = sqrt(r2);
        if( rho < 0.0) h2 = -h2;
        h3 = h1 * h2;
        h7 = exp(-h3 / 2.0);
        if (fabs(rho) < 1.0)
        {
          h6 = fabs(h1 - h2);
          h5 = h6 * h6 / 2.0;
          h6 = h6 / r3;
          AA = 0.5 - h3 / 8.0;
          ab = 3.0 - 2.0 * AA * h5;
          LH = 0.13298076 * h6 * ab * (1 - CND(h6)) - exp(-h5 / r2) * (ab + AA * r2) * 0.053051647;
          for( i = 0 ; i<= 4;i++)
          {
            r1 = r3 * X[i];
            rr = r1 * r1;
            r2 = sqrt(1 - rr);
            if( h7 == 0)
            {
              h8 = 0;
            }
            else
            {
              h8 = exp(-h3 / (1.0 + r2)) / r2 / h7;
            }
            LH = LH - W[i] * exp(-h5 / rr) * (h8 - 1 - AA * rr);
          }
        }

        result = LH * r3 * h7 + CND(Min(h1, h2));
        if (rho < 0.0)
        {
          result = CND(h1) - result;
        }
    }
    else
    {
        h3 = h1 * h2;
        if (rho != 0)
        {
          for( i = 0 ; i<= 4; i++)
          {
            r1 = rho * X[i];
            r2 = 1 - r1 * r1;
            LH = LH + W[i] * exp((r1 * h3 - h12) / r2) / sqrt(r2);
          }
        }
        result = CND(h1) * CND(h2) + rho * LH;
    }

    return result;
    
}

double OptionUtil::CBNDGeneral(int TypeFlag, double X, double y, double rho)
{
    switch (TypeFlag)
        {
            case 1:
                return CBND2(X, y, rho);
            case 2:
                return CBND3(X, y, rho);
            case 3:
                return CBND4(X, y, rho);
            case 4:
                return CBND(X, y, rho);
            default:
                return CBND(X, y, rho);
        }
    
}

/*
 * //http://home.online.no/~pjacklam/notes/invnorm/impl/sprouse/ltqnorm.c
 */
double OptionUtil::CNDEV(double p)
{
    double a[6] =
    {
        -3.969683028665376e+01,
         2.209460984245205e+02,
        -2.759285104469687e+02,
         1.383577518672690e+02,
        -3.066479806614716e+01,
         2.506628277459239e+00
    };

    double b[5] =
    {
        -5.447609879822406e+01,
         1.615858368580409e+02,
        -1.556989798598866e+02,
         6.680131188771972e+01,
        -1.328068155288572e+01
    };

    double c[6] =
    {
        -7.784894002430293e-03,
        -3.223964580411365e-01,
        -2.400758277161838e+00,
        -2.549732539343734e+00,
         4.374664141464968e+00,
         2.938163982698783e+00
    };

    double d[4] =
    {
        7.784695709041462e-03,
        3.224671290700398e-01,
        2.445134137142996e+00,
        3.754408661907416e+00
    };

    double LOW = 0.02425;
    double HIGH = 0.97575;

    double q = 0.0, r = 0.0;
    if (p < 0 || p > 1)
    {
        return 0.0;
    }
    else if (p == 0)
    {
        return -999999.99;
    }
    else if (p == 1)
    {
        return 999999.99;
    }
    else if (p < LOW)
    {
        /* Rational approximation for lower region */
        q = sqrt(-2.0 * log(p));
        return (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) /
            ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0);

    }
    else if (p > HIGH)
    {
        /* Rational approximation for upper region */
        q = sqrt(-2.0 * log(1.0 - p));
        return -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) /
            ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0);

    }
    else
    {
        /* Rational approximation for central region */
        q = p - 0.5;
        r = q * q;
        return (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q /
            (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1.0);


    }


}

double OptionUtil::GBlackScholes(string CallPutFlag, double S, double x,
            double T, double r, double b, double v)
{
    double d1,d2;
    d1 = (log(S / x) + (b + 0.5*v*v) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);

    if(CallPutFlag =="c")
    {
        return S * exp((b - r) * T) * CND(d1) - x * exp(-r * T) * CND(d2);
    }
    else if (CallPutFlag=="p")
    {
        return x * exp(-r * T) * CND(-d2) - S * exp((b - r) * T) * CND(-d1);
    }
    else
    {
        return 0.0;
    }

}

 double OptionUtil::Black76(string CallPutFlag, double F, double x, double T, double r, double v)
{
    double d1,d2;
    d1 = (log(F / x) + (0.5*v*v) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    if(CallPutFlag =="c")
    {
        return exp(-r * T) * (F * CND(d1) - x * CND(d2));
    }
    else if (CallPutFlag =="p")
    {
        return exp(-r * T) * (x * CND(-d2) - F * CND(-d1));
    }
    else
    {
        return 0.0;
    }

}

double OptionUtil::GarmanKolhagen(string CallPutFlag, double S, double x, double T, double r, double rf, double v)
{
    double d1, d2;

    d1 = (log(S / x) + (r - rf + 0.5*v*v) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    if(CallPutFlag =="c")
    {
        return S * exp(-rf * T) * CND(d1) - x * exp(-r * T) * CND(d2);
    }
    else if(CallPutFlag =="p")
    {
        return x * exp(-r * T) * CND(-d2) - S * exp(-rf * T) * CND(-d1);
    }
    else
    {
        return 0.0;
    }
}

double OptionUtil::GDdeltaDvol(double S, double x, double T, double r, double b, double v)
{
    double d1, d2;

    d1 = (log(S / x) + (b + 0.5*v * v ) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    return -exp((b - r) * T) * d2 / v * ND(d1);

}

double OptionUtil::GDdeltaDvolDvol(double S, double x, double T, double r, double b, double v)
{
    double d1,d2;

    d1 = (log(S / x) + (b + 0.5*v * v) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    return GDdeltaDvol(S, x, T, r, b, v) * (1.0 / v) * (d1 * d2 - d1 / d2 - 1);
}

double OptionUtil::GVegaFromDelta(double S, double T, double r, double b, double delta)
{
    return S * exp((b - r) * T) * sqrt(T) * ND(CNDEV(exp((r - b) * T) * fabs(delta)));
}

double OptionUtil::GGammaFromDelta(double S, double T, double r, double b, double v, double delta)
{
    return exp((b - r) * T) * ND(CNDEV(exp((r - b) * T) * fabs(delta))) / (S * v * sqrt(T));
}

double OptionUtil::GRNDFromInTheMoneyProb(double x, double T, double r, double v, double Probability)
{
    return exp(-r * T) * ND(CNDEV(Probability)) / (x * v * sqrt(T));
}

double OptionUtil::GGammaPFromDelta(double S, double T, double r, double b, double v, double delta)
{
    return S *0.01 * GGammaFromDelta(S, T, r, b, v, delta);
}

double OptionUtil::GVegaPFromDelta(double S, double T, double r, double b, double v, double delta)
{
    return v * 0.1 * GVegaFromDelta(S, T, r, b, delta);
}

double OptionUtil::MaxDdeltaDvolAsset(string UpperLowerFlag, double x, double T, double b, double v)
{
  if(UpperLowerFlag == "l")
  {
          return x * exp(-b * T - v * sqrt(T) * sqrt(4.0 + T * v * v) / 2.0);
  }
  else if(UpperLowerFlag=="u")
  {
          return x * exp(-b * T + v * sqrt(T) * sqrt(4.0 + T * v * v) / 2.0);
  }
  else
  {
          return 0.0;
  }

}

double OptionUtil::MaxDdeltaDvolStrike(string UpperLowerFlag, double S, double T, double b, double v)
{
  if(UpperLowerFlag == "l")
  {
          return S * exp(b * T - v * sqrt(T) * sqrt(4.0 + T * v * v) / 2.0);
  }
  else if(UpperLowerFlag == "u")
  {
          return S * exp(b * T + v * sqrt(T) * sqrt(4.0 + T * v * v) / 2.0);
  }
  else
  {
          return 0.0;
  }

}

double OptionUtil::GMaxGammaVegaatX(double S, double b, double T, double v)
{
    return S * exp((b + v * v *0.5) * T);
}

double OptionUtil::GMaxGammaatS(double x, double b, double T, double v)
{
    return x * exp((-b - 3.0 * v * v * 0.5) * T);
}

double OptionUtil::GMaxVegaatS(double x, double b, double T, double v)
{
    return x * exp((-b + v * v * 0.5) * T);
}

double OptionUtil::GForwardDelta(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    double d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    if(CallPutFlag =="c")
    {
        return exp(-r * T) * CND(d1);
    }
    else if(CallPutFlag == "p")
    {
        return exp(-r * T) * (CND(d1) - 1.0);
    }
    else
    {
        return 0.0;
    }

}

double OptionUtil::GDzetaDvol(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    double d1, d2;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    if(CallPutFlag == "c")
    {
        return -ND(d2) * d1 / v;
    }
    else
    {
        return ND(d2) * d1 / v;
    }

}

double OptionUtil::GDzetaDtime(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    double d1, d2;
    d1 = (log(S / x) + (b + v * v *0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    if(CallPutFlag == "c")
    {
        return ND(d2) * (b / (v * sqrt(T)) - d1 / (2.0 * T));
    }
    else
    {
        return -ND(d2) * (b / (v * sqrt(T)) - d1 / (2.0 * T));
    }

}

double OptionUtil::GInTheMoneyProbability(string CallPutFlag, double S, double x, double T, double b, double v)
{
    double d2;
    d2 = (log(S / x) + (b - v * v *0.5) * T) / (v * sqrt(T));
    if(CallPutFlag == "c")
    {
        return  CND(d2);
    }
    else if(CallPutFlag == "p")
    {
        return  CND(-d2);
    }
    else
    {
        return CND(d2);
    }

}

double OptionUtil::GBreakEvenProbability(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
   double d2;
   double xNew;
   if(CallPutFlag == "c")
   {
       xNew = x + GBlackScholes("c", S, x, T, r, b, v) * exp(r * T);
       d2 = (log(S / xNew) + (b - v * v * 0.5) * T) / (v * sqrt(T));
       return CND(d2);
   }
   else if(CallPutFlag == "p")
   {
           xNew = x - GBlackScholes("p", S, x, T, r, b, v) * exp(r * T);
           d2 = (log(S / x) + (b - v * v * 0.5) * T) / (v * sqrt(T));
           return CND(-d2);
   }
   else
   {
           return 0.0;
   }

}

double OptionUtil::GStrikeFromInTheMoneyProb(string CallPutFlag, double S, double v, double T, double b, double InTheMoneyProb)
{
  if(CallPutFlag =="c")
  {
          return  S * exp(-CNDEV(InTheMoneyProb) * v * sqrt(T) + (b - v * v * 0.5) * T);
  }
  else
  {
          return  S * exp(CNDEV(InTheMoneyProb) * v * sqrt(T) + (b - v * v * 0.5) * T);
  }

}

double OptionUtil::GStrikeFromDelta(string CallPutFlag, double S, double T, double r, double b, double v, double delta)
{
  if(CallPutFlag=="c")
  {
          return  S * exp(-CNDEV(delta * exp((r - b) * T)) * v * sqrt(T) + (b + v * v * 0.5) * T);
  }
  else
  {
          return  S * exp(CNDEV(-delta * exp((r - b) * T)) * v * sqrt(T) + (b + v * v * 0.5) * T);
  }

}

double OptionUtil::InTheMoneyProbFromDelta(string CallPutFlag, double S, double T, double r, double b, double v, double delta)
{
  if(CallPutFlag == "c")
  {
          return  CND(CNDEV(delta / exp((b - r) * T)) - v * sqrt(T));
  }
  else
  {
          return  CND(CNDEV(-delta / exp((b - r) * T)) + v * sqrt(T));
  }

}

double OptionUtil::GDeltaFromInTheMoneyProb(string CallPutFlag, double S, double T, double r, double b, double v, double InTheMoneyProb)
{
  if(CallPutFlag == "c")
  {
          return  CND(CNDEV(InTheMoneyProb * exp((b - r) * T)) - v * sqrt(T));
  }
  else
  {
          return  -CND(CNDEV(InTheMoneyProb * exp((b - r) * T)) + v * sqrt(T));
  }

}

double OptionUtil::GDeltaMirrorStrike(double S, double T, double b, double v)
{
  return S * exp((b + v * v * 0.5) * T);
}

double OptionUtil::GProbabilityMirrorStrike(double S, double T, double b, double v)
{
    return S * exp((b - v * v * 0.5) * T);
}

double OptionUtil::GDeltaMirrorCallPutStrike(double S, double x, double T, double b, double v)
{
    return S * S / x * exp((2.0 * b + v * v) * T);
}

double OptionUtil::GGamma(double S, double x, double T, double r, double b, double v)
{
    double d1;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    return exp((b - r) * T) * ND(d1) / (S * v * sqrt(T));
}

double OptionUtil::GSaddleGamma(double x, double T, double r, double b, double v)
{
    return sqrt(exp(1.0) / 3.141592653590) * sqrt((2.0 * b - r) / (v * v) + 1.0) / x;
}

double OptionUtil::GGammaP(double S, double x, double T, double r, double b, double v)
{
    return S * GGamma(S, x, T, r, b, v) * 0.01;
}

double OptionUtil::GDelta(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    double d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    if(CallPutFlag=="c")
    {
        return  exp((b - r) * T) * CND(d1);
    }
    else
    {
        return  -exp((b - r) * T) * CND(-d1);
    }
}

double OptionUtil::GStrikeDelta(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    double d2=(log(S / x) + (b - v * v * 0.5) * T) / (v * sqrt(T));
    if(CallPutFlag =="c")
    {
        return  -exp(-r * T) * CND(d2);
    }
    else
    {
        return  exp(-r * T) * CND(-d2);
    }

}

double OptionUtil::GElasticity(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    return GDelta(CallPutFlag, S, x, T, r, b, v) * S / GBlackScholes(CallPutFlag, S, x, T, r, b, v);
}

double OptionUtil::GDgammaDvol(double S, double x, double T, double r, double b, double v)
{
    double d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    double d2 = d1 - v * sqrt(T);
    return GGamma(S, x, T, r, b, v) * ((d1 * d2 - 1.0) / v);
}

double OptionUtil::GDgammaPDvol(double S, double x, double T, double r, double b, double v)
{
    double d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    double d2 = d1 - v * sqrt(T);
    return S * 0.01 * GGamma(S, x, T, r, b, v) * ((d1 * d2 - 1.0) / v);
}

double OptionUtil::GDgammaDspot(double S, double x, double T, double r, double b, double v)
{
    double d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    return -GGamma(S, x, T, r, b, v) * (1.0 + d1 / (v * sqrt(T))) / S;
}

double OptionUtil::GDgammaPDspot(double S, double x, double T, double r, double b, double v)
{
    double d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    return -GGamma(S, x, T, r, b, v) * (d1) / (100.0 * v * sqrt(T));
}

double OptionUtil::GRiskNeutralDensity(double S, double x, double T, double r, double b, double v)
{
    double d2 = (log(S / x) + (b - v * v * 0.5) * T) / (v * sqrt(T));
    return exp(-r * T) * ND(d2) / (x * v * sqrt(T));
}

double OptionUtil::GTheta(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    double d1,d2;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    if(CallPutFlag =="c")
    {
        return  -S * exp((b - r) * T) * ND(d1) * v / (2.0 * sqrt(T)) - (b - r) * S * exp((b - r) * T) * CND(d1) - r * x * exp(-r * T) * CND(d2);
    }
    else if(CallPutFlag =="p")
    {
        return -S * exp((b - r) * T) * ND(d1) * v / (2.0 * sqrt(T)) + (b - r) * S * exp((b - r) * T) * CND(-d1) + r * x * exp(-r * T) * CND(-d2);
    }
    else
    {
        return 0.0;
    }

}

double OptionUtil::GThetaDriftLess(double S, double x, double T, double r, double b, double v)
{
    double d1;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    return -S * exp((b - r) * T) * ND(d1) * v / (2 * sqrt(T));
}

double OptionUtil::GVarianceVega(double S, double x, double T, double r, double b, double v)
{
    double d1;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    return S * exp((b - r) * T) * ND(d1) * sqrt(T) / (2.0 * v);
}

double OptionUtil::GVarianceVomma(double S, double x, double T, double r, double b, double v)
{
    double d1, d2;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    return S * exp((b - r) * T) * sqrt(T) / (4.0 * v * v * v) * ND(d1) * (d1 * d2 - 1.0);

}

double OptionUtil::GVarianceUltima(double S, double x, double T, double r, double b, double v)
{
    double d1, d2;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    return S * exp((b - r) * T) * sqrt(T) / (8 * pow(v,5)) * ND(d1) * ((d1 * d2 - 1) * (d1 * d2 - 3) - (d1 * d1 + d2 * d2));
}


double OptionUtil::GVarianceDelta(double S, double x, double T, double r, double b, double v)
{
    double d1, d2;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    return S * exp((b - r) * T) * ND(d1) * (-d2) / (2.0 * v * v );

}

double OptionUtil::GVega(double S, double x, double T, double r, double b, double v)
{
    double d1;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    return S * exp((b - r) * T) * ND(d1) * sqrt(T);
}

double OptionUtil::GVegaP(double S, double x, double T, double r, double b, double v)
{
    return v * 0.1 * GVega(S, x, T, r, b, v);
}

double OptionUtil::GDdeltaDtime(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    double d1, d2;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);

    if(CallPutFlag == "c")
    {
        return -exp((b - r) * T) * (ND(d1) * (b / (v * sqrt(T)) - d2 / (2.0 * T)) + (b - r) * CND(d1));
    }
    else if(CallPutFlag =="p")
    {
        return -exp((b - r) * T) * (ND(d1) * (b / (v * sqrt(T)) - d2 / (2.0 * T)) - (b - r) * CND(-d1));
    }
    else
    {
        return 0.0;
    }

}

double OptionUtil::GProfitLossSTD(string TypeFlag, string CallPutFlag, double S, double x, double T, double r, double b, double v, int NHedges)
{
  if(TypeFlag=="a")
  {
          return  sqrt(3.141592653590 / 4.0) * GVega(S, x, T, r, b, v) * v / sqrt(NHedges);
  }
  else if(TypeFlag=="p")
  {
          return  sqrt(3.141592653590 / 4.0) * GVega(S, x, T, r, b, v) * v / sqrt(NHedges) / GBlackScholes(CallPutFlag, S, x, T, r, b, v);
  }
  else
  {
          return 0.0;
  }

}

double OptionUtil::GDvegaDvol(double S, double x, double T, double r, double b, double v)
{
    double d1, d2;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    return GVega(S, x, T, r, b, v) * d1 * d2 / v;
}

double OptionUtil::GDvegaPDvol(double S, double x, double T, double r, double b, double v)
{
    double d1, d2;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    return GVegaP(S, x, T, r, b, v) * d1 * d2 / v;
}

double OptionUtil::GDvegaDtime(double S, double x, double T, double r, double b, double v)
{
    double d1, d2;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    return GVega(S, x, T, r, b, v) * (r - b + b * d1 / (v * sqrt(T)) - (1.0 + d1 * d2) / (2.0 * T));
}

double OptionUtil::GDvommaDvol(double S, double x, double T, double r, double b, double v)
{
    double d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    double d2 = d1 - v * sqrt(T);
    return GDvegaDvol(S, x, T, r, b, v) * (1.0 / v) * (d1 * d2 - d1 / d2 - d2 / d1 - 1.0);

}

double OptionUtil::GDgammaDtime(double S, double x, double T, double r, double b, double v)
{
    double d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    double d2 = d1 - v * sqrt(T);
    return GGamma(S, x, T, r, b, v) * (r - b + b * d1 / (v * sqrt(T)) + (1.0 - d1 * d2) / (2.0 * T));
}

double OptionUtil::GDgammaPDtime(double S, double x, double T, double r, double b, double v)
{
    double d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    double d2 = d1 - v * sqrt(T);
    return GGammaP(S, x, T, r, b, v) * (r - b + b * d1 / (v * sqrt(T)) + (1.0 - d1 * d2) / (2.0 * T));
}

double OptionUtil::GVegaLeverage(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    return GVega(S, x, T, r, b, v) * v / GBlackScholes(CallPutFlag, S, x, T, r, b, v);
}

double OptionUtil::GRho(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    double d1, d2;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    d2 = d1 - v * sqrt(T);
    if(CallPutFlag == "c")
    {
        return T * x * exp(-r * T) * CND(d2);
    }
    else if(CallPutFlag =="p")
    {
        return  -T * x * exp(-r * T) * CND(-d2);
    }
    else
    {
        return 0.0;
    }

}

double OptionUtil::GRhoFO(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    return -T * GBlackScholes(CallPutFlag, S, x, T, r, 0.0, v);
}

double OptionUtil::GPhi(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    double d1;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    if(CallPutFlag=="c")
    {
        return  -T * S * exp((b - r) * T) * CND(d1);
    }
    else if(CallPutFlag=="p")
    {
        return  T * S * exp((b - r) * T) * CND(-d1);
    }
    else
    {
        return 0.0;
    }

}

double OptionUtil::GCarry(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
    double d1;
    d1 = (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    if(CallPutFlag =="c")
    {
        return T * S * exp((b - r) * T) * CND(d1);
    }
    else if(CallPutFlag=="p")
    {
        return -T * S * exp((b - r) * T) * CND(-d1);
    }
    else
    {
        return 0.0;
    }
}

double OptionUtil::GBlackScholesVariance(string CallPutFlag, double S, double x, double T, double r, double b, double v)
{
  double d1, d2;
  d1 = (log(S / x) + (b + v * 0.5) * T) / sqrt(v * T);
  d2 = d1 - sqrt(v * T);
  if(CallPutFlag=="c")
  {
          return S * exp((b - r) * T) * CND(d1) - x * exp(-r * T) * CND(d2);
  }
  else if(CallPutFlag=="p")
  {
          return x * exp(-r * T) * CND(-d2) - S * exp((b - r) * T) * CND(-d1);
  }
  else
  {
          return 0.0;
  }

}

double OptionUtil::GBlackScholesNGreeks(string OutPutFlag, string CallPutFlag, double S, double x, double T, double r, double b, double v, double dS)
{
      if(OutPutFlag==  "p") // Then ' Value
      {
                return  GBlackScholes(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag==  "d")// Then 'Delta
      {
                return  (GBlackScholes(CallPutFlag, S + dS, x, T, r, b, v) - GBlackScholes(CallPutFlag, S - dS, x, T, r, b, v)) / (2.0 * dS);
      }
      else if(OutPutFlag==  "dP")// Then 'Delta
      {
                return  (GBlackScholes(CallPutFlag, S * (1 + dS), x, T, r, b, v) - GBlackScholes(CallPutFlag, S * (1 - dS), x, T, r, b, v)) * 2.0 / S;
      }
      else if(OutPutFlag==  "e")// Then 'Elasticity
      {
                return  (GBlackScholes(CallPutFlag, S + dS, x, T, r, b, v) - GBlackScholes(CallPutFlag, S - dS, x, T, r, b, v)) / (2.0 * dS) * S / GBlackScholes(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag==  "g")// Then 'Gamma
      {
                return (GBlackScholes(CallPutFlag, S + dS, x, T, r, b, v) - 2.0 * GBlackScholes(CallPutFlag, S, x, T, r, b, v) + GBlackScholes(CallPutFlag, S - dS, x, T, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag==  "gv")// Then 'DGammaDVol
      {
                return (GBlackScholes(CallPutFlag, S + dS, x, T, r, b, v + 0.01) - 2.0 * GBlackScholes(CallPutFlag, S, x, T, r, b, v + 0.01) + GBlackScholes(CallPutFlag, S - dS, x, T, r, b, v + 0.01)
              - GBlackScholes(CallPutFlag, S + dS, x, T, r, b, v - 0.01) + 2.0 * GBlackScholes(CallPutFlag, S, x, T, r, b, v - 0.01) - GBlackScholes(CallPutFlag, S - dS, x, T, r, b, v - 0.01)) / (2.0 * 0.01 * dS * dS) / 100.0;
      }
      else if(OutPutFlag==  "gp")// Then 'GammaP
      {
                return  S / 100.0 * (GBlackScholes(CallPutFlag, S + dS, x, T, r, b, v) - 2.0 * GBlackScholes(CallPutFlag, S, x, T, r, b, v) + GBlackScholes(CallPutFlag, S - dS, x, T, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag==  "dddv")//Then 'DDeltaDvol
      {
                return  1.0 / (4.0 * dS * 0.01) * (GBlackScholes(CallPutFlag, S + dS, x, T, r, b, v + 0.01) - GBlackScholes(CallPutFlag, S + dS, x, T, r, b, v - 0.01)
                   - GBlackScholes(CallPutFlag, S - dS, x, T, r, b, v + 0.01) + GBlackScholes(CallPutFlag, S - dS, x, T, r, b, v - 0.01)) / 100.0;
      }
      else if(OutPutFlag==  "v")// Then 'Vega
      {
                return  (GBlackScholes(CallPutFlag, S, x, T, r, b, v + 0.01) - GBlackScholes(CallPutFlag, S, x, T, r, b, v - 0.01)) / 2.0;
      }
      else if(OutPutFlag==  "vp")// Then 'VegaP
      {
                return  v / 0.1 * (GBlackScholes(CallPutFlag, S, x, T, r, b, v + 0.01) - GBlackScholes(CallPutFlag, S, x, T, r, b, v - 0.01)) / 2.0;
      }
      else if(OutPutFlag==  "dvdv")// Then 'DvegaDvol
      {
                return (GBlackScholes(CallPutFlag, S, x, T, r, b, v + 0.01) - 2.0 * GBlackScholes(CallPutFlag, S, x, T, r, b, v) + GBlackScholes(CallPutFlag, S, x, T, r, b, v - 0.01));
      }
      else if(OutPutFlag==  "t")// Then 'Theta
      {
                 if (T <= 1.0 / 256.0)
                 {
                     return GBlackScholes(CallPutFlag, S, x, 0.00001, r, b, v) - GBlackScholes(CallPutFlag, S, x, T, r, b, v);
                 }
                 else
                 {
                     return  GBlackScholes(CallPutFlag, S, x, T - 1.0 / 256.0, r, b, v) - GBlackScholes(CallPutFlag, S, x, T, r, b, v);
                 }
      }
      else if(OutPutFlag==  "r")// Then 'Rho
      {
                return (GBlackScholes(CallPutFlag, S, x, T, r + 0.01, b + 0.01, v) - GBlackScholes(CallPutFlag, S, x, T, r - 0.01, b - 0.01, v)) / (2.0);
      }
      else if(OutPutFlag==  "fr")// Then 'Futures options rho
      {
                return (GBlackScholes(CallPutFlag, S, x, T, r + 0.01, 0.0, v) - GBlackScholes(CallPutFlag, S, x, T, r - 0.01, 0.0, v)) / (2.0);
      }
      else if(OutPutFlag== "f")// Then 'Rho2
      {
                return (GBlackScholes(CallPutFlag, S, x, T, r, b - 0.01, v) - GBlackScholes(CallPutFlag, S, x, T, r, b + 0.01, v)) / (2.0);
      }
      else if(OutPutFlag==  "b")// Then 'Carry
      {
                return (GBlackScholes(CallPutFlag, S, x, T, r, b + 0.01, v) - GBlackScholes(CallPutFlag, S, x, T, r, b - 0.01, v)) / (2.0);
      }
      else if(OutPutFlag==  "s")// Then 'Speed
      {
                return  1.0 / (dS * dS * dS) * (GBlackScholes(CallPutFlag, S + 2.0 * dS, x, T, r, b, v) - 3.0 * GBlackScholes(CallPutFlag, S + dS, x, T, r, b, v)
                                + 3.0 * GBlackScholes(CallPutFlag, S, x, T, r, b, v) - GBlackScholes(CallPutFlag, S - dS, x, T, r, b, v));
      }
      else if(OutPutFlag==  "dx")// Then 'Strike Delta
      {
                return (GBlackScholes(CallPutFlag, S, x + dS, T, r, b, v) - GBlackScholes(CallPutFlag, S, x - dS, T, r, b, v)) / (2.0 * dS);
      }
      else if(OutPutFlag== "dxdx")//Then 'Gamma
      {
                return (GBlackScholes(CallPutFlag, S, x + dS, T, r, b, v) - 2.0 * GBlackScholes(CallPutFlag, S, x, T, r, b, v) + GBlackScholes(CallPutFlag, S, x - dS, T, r, b, v)) / (dS * dS);
      }
      else
      {
              return 0.0;
      }
}

double OptionUtil::CGBlackScholes(string OutPutFlag, string CallPutFlag, double S, double x, double T, double r, double b, double v, double delta, double InTheMoneyProb, double ThetaDays)
{
      if(OutPutFlag==  "p")// Then 'Value
      {
                return  GBlackScholes(CallPutFlag, S, x, T, r, b, v);
      }
      /////DELTA GREEKS
      else if(OutPutFlag==  "d")// Then 'Delta
      {
                return  GDelta(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag== "df")// Then 'Forward Delta
      {
                return GForwardDelta(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag== "dddv")// Then 'DDeltaDvol
      {
                return GDdeltaDvol(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag==  "dvv")// Then 'DDeltaDvolDvol
      {
                return GDdeltaDvolDvol(S, x, T, r, b, v) / 10000.0;
      }
      else if(OutPutFlag==  "dt")// Then 'DDeltaDtime/Charm
      {
                return GDdeltaDtime(CallPutFlag, S, x, T, r, b, v) / 256.0;
      }
      else if(OutPutFlag==  "dmx")// Then
      {
                return S * S / x * exp((2.0 * b + v * v) * T);
      }
      else if(OutPutFlag==  "e")// Then ' Elasticity
      {
                 return GElasticity(CallPutFlag, S, x, T, r, b, v);
      }
      //////GAMMA GREEKS
      else if(OutPutFlag==  "sg")// Then 'SaddleGamma
      {
                return  GSaddleGamma(x, T, r, b, v);
      }
      else if(OutPutFlag==  "g")// Then 'Gamma
      {
                return  GGamma(S, x, T, r, b, v);
      }
      else if(OutPutFlag==  "s")// Then 'DgammaDspot/speed
      {
                return  GDgammaDspot(S, x, T, r, b, v);
      }
      else if(OutPutFlag==  "gv")// Then 'DgammaDvol/Zomma
      {
                return  GDgammaDvol(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag==  "gt")// Then 'DgammaDtime
      {
                return  GDgammaDtime(S, x, T, r, b, v) / 256.0;
      }
      else if(OutPutFlag==  "gp")// Then 'GammaP
      {
                return  GGammaP(S, x, T, r, b, v);
      }
      else if(OutPutFlag==  "gps")// Then 'DgammaPDspot
      {
                return GDgammaPDspot(S, x, T, r, b, v);
      }
      else if(OutPutFlag==  "gpv")// Then 'DgammaDvol/Zomma
      {
                return GDgammaPDvol(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag==  "gpt")// Then 'DgammaPDtime
      {
                return GDgammaPDtime(S, x, T, r, b, v) / 256.0;
      }
      //////VEGA GREEKS
      else if(OutPutFlag== "v")// Then 'Vega
      {
                return  GVega(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag==  "vt")// Then 'DvegaDtime
      {
                return GDvegaDtime(S, x, T, r, b, v) / 256.0;
      }
      else if(OutPutFlag==  "dvdv")// Then 'DvegaDvol/Vomma
      {
                return GDvegaDvol(S, x, T, r, b, v) / 10000.0;
      }
      else if(OutPutFlag== "vvv")// Then 'DvommaDvol
      {
                return GDvommaDvol(S, x, T, r, b, v) / 1000000.0;
      }
      else if(OutPutFlag== "vp")// Then 'VegaP
      {
                return GVegaP(S, x, T, r, b, v);
      }
      else if(OutPutFlag== "vpv")// Then 'DvegaPDvol/VommaP
      {
                return GDvegaPDvol(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "vl")// Then 'Vega Leverage
      {
                return GVegaLeverage(CallPutFlag, S, x, T, r, b, v);
      }
      //////VARIANCE GREEKS
      else if(OutPutFlag== "varvega")// Then 'Variance-Vega
      {
                return GVarianceVega(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "vardelta")// Then 'Variance-delta
      {
                return GVarianceDelta(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "varvar")// Then 'Variance-vomma
      {
                return GVarianceVomma(S, x, T, r, b, v) / 10000.0;
      }
      else if(OutPutFlag== "varult")// Then 'Variance-ultima
      {
                return GVarianceUltima(S, x, T, r, b, v) / 1000000.0;
      }
      ///////THETA GREEKS
      else if(OutPutFlag== "t")// Then 'Theta
      {
                return GTheta(CallPutFlag, S, x, T, r, b, v) / 256.0;
      }
      else if(OutPutFlag== "Dlt")// Then 'Drift-less Theta
      {
                return GThetaDriftLess(S, x, T, r, b, v) / 256.0;
      }
      //////RATE/CARRY GREEKS
      else if(OutPutFlag== "r")// Then 'Rho
      {
                return GRho(CallPutFlag, S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "fr")// Then 'Rho futures option
      {
                return GRhoFO(CallPutFlag, S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "b")// Then 'Carry Rho
      {
                return GCarry(CallPutFlag, S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "f")// Then 'Phi/Rho2
      {
                return GPhi(CallPutFlag, S, x, T, r, b, v) / 100.0;
      }
      //////PROB GREEKS
      else if(OutPutFlag== "z")// Then 'Zeta/In-the-money risk neutral probability
      {
                return GInTheMoneyProbability(CallPutFlag, S, x, T, b, v);
      }
      else if(OutPutFlag== "zv")// Then 'DzetaDvol
      {
                return GDzetaDvol(CallPutFlag, S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "zt")// Then 'DzetaDtime
      {
                return GDzetaDtime(CallPutFlag, S, x, T, r, b, v) / 256.0;
      }
      else if(OutPutFlag== "bp")// Then 'Brak even probability
      {
                return GBreakEvenProbability(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag== "dx")// Then 'StrikeDelta
      {
                return GStrikeDelta(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag== "dxdx")// Then 'Risk Neutral Density
      {
                return GRiskNeutralDensity(S, x, T, r, b, v);
      }
      //////FROM DELTA GREEKS
      else if(OutPutFlag== "gfd")// Then 'Gamma from delta
      {
                return GGammaFromDelta(S, T, r, b, v, delta);
      }
      else if(OutPutFlag== "gpfd")// Then 'GammaP from delta
      {
                return GGammaPFromDelta(S, T, r, b, v, delta);
      }
      else if(OutPutFlag== "vfd")// Then 'Vega from delta
      {
                return GVegaFromDelta(S, T, r, b, delta) / 100.0;
      }
      else if(OutPutFlag== "vpfd")// Then 'VegaP from delta
      {
                return GVegaPFromDelta(S, T, r, b, v, delta);
      }
      else if(OutPutFlag== "xfd")// Then 'Strike from delta
      {
                return GStrikeFromDelta(CallPutFlag, S, T, r, b, v, delta);
      }
      else if(OutPutFlag== "ipfd")// Then 'In-the-money probability from delta
      {
                return InTheMoneyProbFromDelta(CallPutFlag, S, T, r, b, v, delta);
      }
      //////FROM IN-THE GREEKS
      else if(OutPutFlag== "xfip")// Then 'Strike from in-the-money probability
      {
                return GStrikeFromInTheMoneyProb(CallPutFlag, S, v, T, b, InTheMoneyProb);
      }
      else if(OutPutFlag== "RNDfip")// Then 'Risk Neutral Density from in-the-money probability
      {
                return GRNDFromInTheMoneyProb(x, T, r, v, InTheMoneyProb);
      }
      else if(OutPutFlag== "dfip")// Then 'Strike from in-the-money probability
      {
                return GDeltaFromInTheMoneyProb(CallPutFlag, S, T, r, b, v, InTheMoneyProb);
      }
      ///////CALCULATIONS
      else if(OutPutFlag== "d1")// Then 'd1
      {
                return (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
      }
      else if(OutPutFlag == "d2")// Then 'd2
      {
                return (log(S / x) + (b - v * v * 0.5) * T) / (v * sqrt(T));
      }
      else if(OutPutFlag== "nd1")// Then 'n(d1)
      {
                return ND((log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T)));
      }
      else if(OutPutFlag== "nd2")// Then 'n(d2)
      {
                return ND((log(S / x) + (b - v * v * 0.5) * T) / (v * sqrt(T)));
      }
      else if(OutPutFlag== "CNDd1")// Then 'N(d1)
      {
                return CND((log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T)));
      }
      else if(OutPutFlag== "CNDd2")// Then 'N(d2)
      {
                return CND((log(S / x) + (b - v * v * 0.5) * T) / (v * sqrt(T)));
      }
      else
      {
                return 0.0;
      }
}

double OptionUtil::EGBlackScholes(string OutPutFlag, string CallPutFlag, double S, double x, double T, double r, double b, double v, double delta, double InTheMoneyProb, double ThetaDays)
{
      if(OutPutFlag == "p")// Then 'Value
      {
                return GBlackScholes(CallPutFlag, S, x, T, r, b, v);
      }
      ////Delta greeks
      else if(OutPutFlag== "d")// Then 'Delta
      {
                return GDelta(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag== "df")// Then 'Forward Delta
      {
                return GForwardDelta(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag== "dddv")// Then 'DDeltaDvol
      {
                return GDdeltaDvol(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag == "dvv")// Then 'DDeltaDvolDvol
      {
                return GDdeltaDvolDvol(S, x, T, r, b, v) / 10000.0;
      }
      else if(OutPutFlag == "dt")// Then 'DDeltaDtime/Charm
      {
                return GDdeltaDtime(CallPutFlag, S, x, T, r, b, v) / 256.0;
      }
      else if(OutPutFlag == "dmx")// Then
      {
                return S * S / x * exp((2.0 * b + v * v) * T);
      }
      else if(OutPutFlag == "e")// Then ' Elasticity
      {
                return GElasticity(CallPutFlag, S, x, T, r, b, v);
      }
      /////GAMMA GREEKS
      else if(OutPutFlag == "sg")// Then 'SaddleGamma
      {
                return GSaddleGamma(x, T, r, b, v);
      }
      else if(OutPutFlag == "g")// Then 'Gamma
      {
                return GGamma(S, x, T, r, b, v);
      }
      else if(OutPutFlag == "s")// Then 'DgammaDspot/speed
      {
                return GDgammaDspot(S, x, T, r, b, v);
      }
      else if(OutPutFlag== "gv")// Then 'DgammaDvol/Zomma
      {
                return GDgammaDvol(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "gt")// Then 'DgammaDtime
      {
                return GDgammaDtime(S, x, T, r, b, v) / 256.0;
      }
      else if(OutPutFlag == "gp")// Then 'GammaP
      {
                return GGammaP(S, x, T, r, b, v);
      }
      else if(OutPutFlag == "gps")// Then 'DgammaPDspot
      {
                return GDgammaPDspot(S, x, T, r, b, v);
      }
      else if(OutPutFlag == "gpv")// Then 'DgammaDvol/Zomma
      {
                return GDgammaPDvol(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "gpt")// Then 'DgammaPDtime
      {
                return GDgammaPDtime(S, x, T, r, b, v) / 256.0;
      }
      /////VEGA GREEKS
      else if(OutPutFlag == "v")// Then 'Vega
      {
                return GVega(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag == "vt")// Then 'DvegaDtime
      {
                return GDvegaDtime(S, x, T, r, b, v) / 256.0;
      }
      else if(OutPutFlag== "dvdv")// Then 'DvegaDvol/Vomma
      {
                return GDvegaDvol(S, x, T, r, b, v) / 10000.0;
      }
      else if(OutPutFlag == "vvv")// Then 'DvommaDvol
      {
                return GDvommaDvol(S, x, T, r, b, v) / 1000000.0;
      }
      else if(OutPutFlag == "vp")// Then 'VegaP
      {
                return GVegaP(S, x, T, r, b, v);
      }
      else if(OutPutFlag == "vpv")// Then 'DvegaPDvol/VommaP
      {
                return GDvegaPDvol(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag == "vl")// Then 'Vega Leverage
      {
                return GVegaLeverage(CallPutFlag, S, x, T, r, b, v);
      }
      //////VARIANCE GREEKS
      else if(OutPutFlag == "varvega")// Then 'Variance-Vega
      {
                return GVarianceVega(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag == "vardelta")// Then 'Variance-delta
      {
                return GVarianceDelta(S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag == "varvar")// Then 'Variance-vomma
      {
                return GVarianceVomma(S, x, T, r, b, v) / 10000.0;
      }
      else if(OutPutFlag == "varult")// Then 'Variance-ultima
      {
                return GVarianceUltima(S, x, T, r, b, v) / 1000000.0;
      }
      //////THETA GREEKS
      else if(OutPutFlag== "t")// Then 'Theta
      {
                return GTheta(CallPutFlag, S, x, T, r, b, v) / 256.0;
      }
      else if(OutPutFlag== "Dlt")// Then 'Drift-less Theta
      {
                return GThetaDriftLess(S, x, T, r, b, v) / 256.0;
      }
      //////RATE/CARRY GREEKS
      else if(OutPutFlag == "r")// Then 'Rho
      {
                return GRho(CallPutFlag, S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "fr")// Then 'Rho futures option
      {
                return GRhoFO(CallPutFlag, S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "b")// Then 'Carry Rho
      {
                return GCarry(CallPutFlag, S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "f")// Then 'Phi/Rho2
      {
                return GPhi(CallPutFlag, S, x, T, r, b, v) / 100.0;
      }
      //////PROB GREEKS
      else if(OutPutFlag== "z")// Then 'Zeta/In-the-money risk neutral probability
      {
                return GInTheMoneyProbability(CallPutFlag, S, x, T, b, v);
      }
      else if(OutPutFlag== "zv")// Then 'DzetaDvol
      {
                return GDzetaDvol(CallPutFlag, S, x, T, r, b, v) / 100.0;
      }
      else if(OutPutFlag== "zt")// Then 'DzetaDtime
      {
                return GDzetaDtime(CallPutFlag, S, x, T, r, b, v) / 256.0;
      }
      else if(OutPutFlag== "bp")// Then 'Brak even probability
      {
                return GBreakEvenProbability(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag== "dx")// Then 'StrikeDelta
      {
                return GStrikeDelta(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag== "dxdx")// Then 'Risk Neutral Density
      {
                return GRiskNeutralDensity(S, x, T, r, b, v);
      }
      //////FROM DELTA GREEKS
      else if(OutPutFlag== "gfd")// Then 'Gamma from delta
      {
                return GGammaFromDelta(S, T, r, b, v, delta);
      }
      else if(OutPutFlag== "gpfd")// Then 'GammaP from delta
      {
                return GGammaPFromDelta(S, T, r, b, v, delta);
      }
      else if(OutPutFlag== "vfd")// Then 'Vega from delta
      {
                return GVegaFromDelta(S, T, r, b, delta) / 100.0;
      }
      else if(OutPutFlag== "vpfd")// Then 'VegaP from delta
      {
                return GVegaPFromDelta(S, T, r, b, v, delta);
      }
      else if(OutPutFlag== "xfd")// Then 'Strike from delta
      {
                return GStrikeFromDelta(CallPutFlag, S, T, r, b, v, delta);
      }
      else if(OutPutFlag== "ipfd")// Then 'In-the-money probability from delta
      {
                return InTheMoneyProbFromDelta(CallPutFlag, S, T, r, b, v, delta);
      }
      //////FROM IN-THE GREEKS
      else if(OutPutFlag== "xfip")// Then 'Strike from in-the-money probability
      {
                return GStrikeFromInTheMoneyProb(CallPutFlag, S, v, T, b, InTheMoneyProb);
      }
      else if(OutPutFlag== "RNDfip")// Then 'Risk Neutral Density from in-the-money probability
      {
                return GRNDFromInTheMoneyProb(x, T, r, v, InTheMoneyProb);
      }
      else if(OutPutFlag== "dfip")// Then 'Strike from in-the-money probability
      {
                return GDeltaFromInTheMoneyProb(CallPutFlag, S, T, r, b, v, InTheMoneyProb);
      }
      //////CALCULATIONS
      else if(OutPutFlag== "d1")// Then 'd1
      {
                return (log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T));
      }
      else if(OutPutFlag== "d2")// Then 'd2
      {
                return (log(S / x) + (b - v * v * 0.5) * T) / (v * sqrt(T));
      }
      else if(OutPutFlag== "nd1")// Then 'n(d1)
      {
                return ND((log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T)));
      }
      else if(OutPutFlag== "nd2")// Then 'n(d2)
      {
                return ND((log(S / x) + (b - v * v * 0.5) * T) / (v * sqrt(T)));
      }
      else if(OutPutFlag== "CNDd1")// Then 'N(d1)
      {
                return CND((log(S / x) + (b + v * v * 0.5) * T) / (v * sqrt(T)));
      }
      else if(OutPutFlag== "CNDd2")// Then 'N(d2)
      {
                return CND((log(S / x) + (b - v * v * 0.5) * T) / (v * sqrt(T)));
      }
      else
      {
              return 0.0;
      }
}

double OptionUtil::GBlackScholesVarianceNGreeks(string OutPutFlag, string CallPutFlag, double S, double x, double T, double r, double b, double v, double dS)
{
      if(OutPutFlag=="p") // Value
      {
              return  GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v);
      }
      else if (OutPutFlag== "d" ) //Then 'Delta
      {
        return  (GBlackScholesVariance(CallPutFlag, S + dS, x, T, r, b, v) - GBlackScholesVariance(CallPutFlag, S - dS, x, T, r, b, v)) / (2.0 * dS);
      }
      else if(OutPutFlag==  "e") //Then 'Elasticity
      {
        return  (GBlackScholesVariance(CallPutFlag, S + dS, x, T, r, b, v) - GBlackScholesVariance(CallPutFlag, S - dS, x, T, r, b, v)) / (2.0 * dS) * S / GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v);
      }
      else if(OutPutFlag =="g") // Then 'Gamma
      {
        return  (GBlackScholesVariance(CallPutFlag, S + dS, x, T, r, b, v) - 2.0 * GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v) + GBlackScholesVariance(CallPutFlag, S - dS, x, T, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag == "gv") // Then 'DGammaDvariance
      {
        return (GBlackScholesVariance(CallPutFlag, S + dS, x, T, r, b, v + 0.01) - 2.0 * GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v + 0.01) + GBlackScholesVariance(CallPutFlag, S - dS, x, T, r, b, v + 0.01)
           - GBlackScholesVariance(CallPutFlag, S + dS, x, T, r, b, v - 0.01) + 2.0 * GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v - 0.01) - GBlackScholesVariance(CallPutFlag, S - dS, x, T, r, b, v - 0.01)) / (2.0 * 0.01 * dS * dS) / 100.0;
      }
      else if(OutPutFlag == "gp") // Then 'GammaP
      {
        return 0.01 * S * (GBlackScholesVariance(CallPutFlag, S + dS, x, T, r, b, v) - 2.0 * GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v) + GBlackScholesVariance(CallPutFlag, S - dS, x, T, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag == "dddv") // Then 'DDeltaDvariance
      {
        return 1.0 / (4.0 * dS * 0.01) * (GBlackScholesVariance(CallPutFlag, S + dS, x, T, r, b, v + 0.01) - GBlackScholesVariance(CallPutFlag, S + dS, x, T, r, b, v - 0.01)
          - GBlackScholesVariance(CallPutFlag, S - dS, x, T, r, b, v + 0.01) + GBlackScholesVariance(CallPutFlag, S - dS, x, T, r, b, v - 0.01)) / 100.0;
      }
      else if(OutPutFlag == "v") // Then 'Variance Vega
      {
         return (GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v + 0.01) - GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v - 0.01)) / 2.0;
      }
      else if(OutPutFlag ==  "vp") // Then 'Variance VegaP
      {
         return v / 0.1 * (GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v + 0.01) - GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v - 0.01)) / 2.0;
      }
      else if(OutPutFlag == "dvdv") //Then 'Variance Dvegavariance
      {
        return (GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v + 0.01) - 2.0 * GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v) + GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v - 0.01));
      }
      else if(OutPutFlag == "t" ) //Then 'Theta
      {
        if (T <= 1.0 / 256.0)
        {
                return GBlackScholesVariance(CallPutFlag, S, x, 0.00001, r, b, v) - GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v);
        }
        else
        {
                return GBlackScholesVariance(CallPutFlag, S, x, T - 1.0 / 256.0, r, b, v) - GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v);
        }
      }
      else if(OutPutFlag== "r") // Then 'Rho
      {
         return (GBlackScholesVariance(CallPutFlag, S, x, T, r + 0.01, b + 0.01, v) - GBlackScholesVariance(CallPutFlag, S, x, T, r - 0.01, b - 0.01, v)) / (2.0);
      }
      else if(OutPutFlag== "fr") // Then 'Futures options rho
      {
         return (GBlackScholesVariance(CallPutFlag, S, x, T, r + 0.01, 0.0, v) - GBlackScholesVariance(CallPutFlag, S, x, T, r - 0.01, 0.0, v)) / (2.0);
      }
      else if(OutPutFlag == "f") // Then 'Rho2
      {
         return (GBlackScholesVariance(CallPutFlag, S, x, T, r, b - 0.01, v) - GBlackScholesVariance(CallPutFlag, S, x, T, r, b + 0.01, v)) / (2.0);
      }
      else if(OutPutFlag == "b") // Then 'Carry
      {
        return (GBlackScholesVariance(CallPutFlag, S, x, T, r, b + 0.01, v) - GBlackScholesVariance(CallPutFlag, S, x, T, r, b - 0.01, v)) / (2.0);
      }
      else if(OutPutFlag== "s") // Then 'Speed
      {
        return 1.0 / (dS * dS * dS) * (GBlackScholesVariance(CallPutFlag, S + 2.0 * dS, x, T, r, b, v) - 3.0 * GBlackScholesVariance(CallPutFlag, S + dS, x, T, r, b, v)
                   + 3.0 * GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v) - GBlackScholesVariance(CallPutFlag, S - dS, x, T, r, b, v));
      }
      else if(OutPutFlag== "dx") // Then 'Strike Delta
      {
         return (GBlackScholesVariance(CallPutFlag, S, x + dS, T, r, b, v) - GBlackScholesVariance(CallPutFlag, S, x - dS, T, r, b, v)) / (2.0 * dS);
      }
      else if(OutPutFlag== "dxdx") // Then 'Gamma
      {
        return (GBlackScholesVariance(CallPutFlag, S, x + dS, T, r, b, v) - 2.0 * GBlackScholesVariance(CallPutFlag, S, x, T, r, b, v) + GBlackScholesVariance(CallPutFlag, S, x - dS, T, r, b, v)) / (dS * dS);
      }
      else
      {
              return 0.0;
      }

}

//implied volatilities for european options
double OptionUtil::GBlackScholesImpVolBisection(string CallPutFlag, double S,
            double x,double T,double r,double b,double cm)
{
    double vLow, vHigh, vi;
    double cLow, cHigh, epsilon;
    int counter;

    vLow = 0.005;
    vHigh = 4;
   // epsilon = 0.00000001;
    epsilon = 0.0001;

    cLow = GBlackScholes(CallPutFlag, S, x, T, r, b, vLow);
    cHigh = GBlackScholes(CallPutFlag, S, x, T, r, b, vHigh);
    counter = 0;
    vi = vLow + (cm - cLow) * (vHigh - vLow) / (cHigh - cLow);
    while (fabs(cm - GBlackScholes(CallPutFlag, S, x, T, r, b, vi)) > epsilon)
    {
        counter = counter + 1;
        if (counter == 500) return 0.0;// return -99999.99 ;
        if (GBlackScholes(CallPutFlag, S, x, T, r, b, vi) < cm)
        {
            vLow = vi;
        }
        else
        {
            vHigh = vi;
        }
        cLow = GBlackScholes(CallPutFlag, S, x, T, r, b, vLow);
        cHigh = GBlackScholes(CallPutFlag, S, x, T, r, b, vHigh);
        vi = vLow + (cm - cLow) * (vHigh - vLow) / (cHigh - cLow);
    }
    return vi;

}

double OptionUtil::GImpliedVolatilityNR(string CallPutFlag, double S, double x,
            double T, double r, double b, double cm, double epsilon)
{
        double vi, ci;
        double vegai;
        double minDiff;

        //Manaster and Koehler seed value (vi)
        vi = sqrt(fabs(log(S / x) + r * T) * 2.0 / T);
        ci = GBlackScholes(CallPutFlag, S, x, T, r, b, vi);
        vegai = GVega(S, x, T, r, b, vi);
        minDiff = fabs(cm - ci);

        while( fabs(cm - ci) >= epsilon && fabs(cm - ci) <= minDiff)
        {
            vi = vi - (ci - cm) / vegai;
            ci = GBlackScholes(CallPutFlag, S, x, T, r, b, vi);
            vegai = GVega(S, x, T, r, b, vi);
            minDiff = fabs(cm - ci);
        }

        if (fabs(cm - ci) < epsilon)
        {
            return vi;
        }
        else{
            return 0.0;
        }
}

//perpetual options
double OptionUtil::PerpetualOption(string CallPutFlag, double S, double X, double r, double b, double v)
{
    double y1, y2;
    //double h;

    y1 = 1.0 / 2.0 - b / (v * v) + sqrt(pow((b / (v * v) - 1.0 / 2.0), 2) + 2.0 * r / (v * v));
    y2 = 1.0 / 2.0 - b / (v * v) - sqrt(pow((b / (v * v) - 1.0 / 2.0), 2) + 2.0 * r / (v * v));
    if(CallPutFlag=="c")
    {
        return X / (y1 - 1.0) * pow(((y1 - 1.0) / y1 * S / X), y1);
    }
    else if(CallPutFlag=="p")
    {
        return X / (1.0 - y2) * pow(((y2 - 1.0) / y2 * S / X), y2);
    }
    else
    {
        return 0.0;
    }

}

double OptionUtil::EPerpetualOption(string OutPutFlag, string CallPutFlag, double S, double X,
            double r, double b, double v, double dS)
{
      if(OutPutFlag== "p")// Then ' Value
      {
                return PerpetualOption(CallPutFlag, S, X, r, b, v);
      }
      else if(OutPutFlag== "d")// Then 'Delta
      {
                return (PerpetualOption(CallPutFlag, S + dS, X, r, b, v) - PerpetualOption(CallPutFlag, S - dS, X, r, b, v)) / (2.0 * dS);
      }
      else if(OutPutFlag== "e")// Then 'Elasticity
      {
                return (PerpetualOption(CallPutFlag, S + dS, X, r, b, v) - PerpetualOption(CallPutFlag, S - dS, X, r, b, v)) / (2.0 * dS) * S / PerpetualOption(CallPutFlag, S, X, r, b, v);
      }
      else if(OutPutFlag== "g")// Then 'Gamma
      {
                return (PerpetualOption(CallPutFlag, S + dS, X, r, b, v) - 2.0 * PerpetualOption(CallPutFlag, S, X, r, b, v) + PerpetualOption(CallPutFlag, S - dS, X, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag== "gv")// Then 'DGammaDVol
      {
                return (PerpetualOption(CallPutFlag, S + dS, X, r, b, v + 0.01) - 2.0 * PerpetualOption(CallPutFlag, S, X, r, b, v + 0.01) + PerpetualOption(CallPutFlag, S - dS, X, r, b, v + 0.01)
                   - PerpetualOption(CallPutFlag, S + dS, X, r, b, v - 0.01) + 2.0 * PerpetualOption(CallPutFlag, S, X, r, b, v - 0.01) - PerpetualOption(CallPutFlag, S - dS, X, r, b, v - 0.01)) / (2.0 * 0.01 * dS * dS) / 100.0;
      }
      else if(OutPutFlag== "gp")// Then 'GammaP
      {
                return S / 100.0 * (PerpetualOption(CallPutFlag, S + dS, X, r, b, v) - 2.0 * PerpetualOption(CallPutFlag, S, X, r, b, v) + PerpetualOption(CallPutFlag, S - dS, X, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag== "dddv")// Then 'DDeltaDvol
      {
                return 1.0 / (4.0 * dS * 0.01) * (PerpetualOption(CallPutFlag, S + dS, X, r, b, v + 0.01) - PerpetualOption(CallPutFlag, S + dS, X, r, b, v - 0.01)
                   - PerpetualOption(CallPutFlag, S - dS, X, r, b, v + 0.01) + PerpetualOption(CallPutFlag, S - dS, X, r, b, v - 0.01)) / 100.0;
      }
      else if(OutPutFlag== "v")// Then 'Vega
      {
                return (PerpetualOption(CallPutFlag, S, X, r, b, v + 0.01) - PerpetualOption(CallPutFlag, S, X, r, b, v - 0.01)) / 2.0;
      }
      else if(OutPutFlag== "vv")// Then 'DvegaDvol/vomma
      {
                return (PerpetualOption(CallPutFlag, S, X, r, b, v + 0.01) - 2.0 * PerpetualOption(CallPutFlag, S, X, r, b, v) + PerpetualOption(CallPutFlag, S, X, r, b, v - 0.01)) / (0.01 * 0.01) / 10000;
      }
      else if(OutPutFlag== "vp")// Then 'VegaP
      {
                return v / 0.1 * (PerpetualOption(CallPutFlag, S, X, r, b, v + 0.01) - PerpetualOption(CallPutFlag, S, X, r, b, v - 0.01)) / 2.0;
      }
      else if(OutPutFlag== "dvdv")// Then 'DvegaDvol
      {
                return (PerpetualOption(CallPutFlag, S, X, r, b, v + 0.01) - 2.0 * PerpetualOption(CallPutFlag, S, X, r, b, v) + PerpetualOption(CallPutFlag, S, X, r, b, v - 0.01));
      }
      else if(OutPutFlag== "r")// Then 'Rho
      {
                return (PerpetualOption(CallPutFlag, S, X, r + 0.01, b + 0.01, v) - PerpetualOption(CallPutFlag, S, X, r - 0.01, b - 0.01, v)) / (2.0);
      }
      else if(OutPutFlag== "fr")// Then 'Futures options rho
      {
                return (PerpetualOption(CallPutFlag, S, X, r + 0.01, b, v) - PerpetualOption(CallPutFlag, S, X, r - 0.01, b, v)) / (2.0);
      }
      else if(OutPutFlag== "f")// Then 'Rho2
      {
                return (PerpetualOption(CallPutFlag, S, X, r, b - 0.01, v) - PerpetualOption(CallPutFlag, S, X, r, b + 0.01, v)) / (2.0);
      }
      else if(OutPutFlag== "b")// Then 'Carry
      {
                return (PerpetualOption(CallPutFlag, S, X, r, b + 0.01, v) - PerpetualOption(CallPutFlag, S, X, r, b - 0.01, v)) / (2.0);
      }
      else if(OutPutFlag== "s")// Then 'Speed
      {
                return 1.0 / (dS * dS * dS) * (PerpetualOption(CallPutFlag, S + 2.0 * dS, X, r, b, v) - 3.0 * PerpetualOption(CallPutFlag, S + dS, X, r, b, v)
                              + 3.0 * PerpetualOption(CallPutFlag, S, X, r, b, v) - PerpetualOption(CallPutFlag, S - dS, X, r, b, v));
      }
      else if(OutPutFlag== "dx")// Then 'Strike Delta
      {
                return (PerpetualOption(CallPutFlag, S, X + dS, r, b, v) - PerpetualOption(CallPutFlag, S, X - dS, r, b, v)) / (2.0 * dS);
      }
      else if(OutPutFlag== "dxdx")// Then 'Strike Gamma
      {
                return (PerpetualOption(CallPutFlag, S, X + dS, r, b, v) - 2.0 * PerpetualOption(CallPutFlag, S, X, r, b, v) + PerpetualOption(CallPutFlag, S, X - dS, r, b, v)) / (dS * dS);
      }
      else
      {
                return 0.0;
      }
}

//american options
double OptionUtil::Kc(double X, double T, double r, double b, double v)
{
    double N, m;
    double su, Si;
    double h2, k;
    double d1, Q2, q2u;
    double LHS, RHS;
    double bi, E;

    // Calculation of seed value, Si
    N = 2.0 * b / (v * v);
    m = 2.0 * r / (v * v);
    q2u = (-(N - 1.0) + sqrt((N - 1.0) * (N - 1.0) + 4.0 * m)) / 2.0;
    su = X / (1.0 - 1.0 / q2u);
    h2 = -(b * T + 2.0 * v * sqrt(T)) * X / (su - X);
    Si = X + (su - X) * (1.0 - exp(h2));

    k = 2.0 * r / (v * v * (1.0 - exp(-r * T)));
    d1 = (log(Si / X) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    Q2 = (-(N - 1.0) + sqrt((N - 1.0) * (N - 1.0) + 4.0 * k)) / 2.0;
    LHS = Si - X;
    RHS = GBlackScholes("c", Si, X, T, r, b, v) + (1.0 - exp((b - r) * T) * CND(d1)) * Si / Q2;
    bi = exp((b - r) * T) * CND(d1) * (1.0 - 1.0 / Q2) + (1.0 - exp((b - r) * T) * CND(d1) / (v * sqrt(T))) / Q2;
    E = 0.000001;
    // Newton Raphson algorithm for finding critical price Si
    while (fabs(LHS - RHS) / X > E)
    {
        Si = (X + RHS - bi * Si) / (1.0 - bi);
        d1 = (log(Si / X) + (b + v * v * 0.5) * T) / (v * sqrt(T));
        LHS = Si - X;
        RHS = GBlackScholes("c", Si, X, T, r, b, v) + (1.0 - exp((b - r) * T) * CND(d1)) * Si / Q2;
        bi = exp((b - r) * T) * CND(d1) * (1.0 - 1.0 / Q2) + (1 - exp((b - r) * T) * ND(d1) / (v * sqrt(T))) / Q2;
    }
    return Si;
}

double OptionUtil::Kp(double X, double T, double r, double b, double v)
{
    double N, m;
    double su, Si;
    double h1, k;
    double d1, q1u, Q1;
    double LHS, RHS;
    double bi, E;

    //    '// Calculation of seed value, Si
    N = 2 * b / (v * v);
    m = 2 * r / (v * v);
    q1u = (-(N - 1.0) - sqrt((N - 1.0) * (N - 1.0) + 4.0 * m)) / 2.0;
    su = X / (1 - 1 / q1u);
    h1 = (b * T - 2 * v * sqrt(T)) * X / (X - su);
    Si = su + (X - su) * exp(h1);


    k = 2.0 * r / (v * v * (1.0 - exp(-r * T)));
    d1 = (log(Si / X) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    Q1 = (-(N - 1.0) - sqrt((N - 1.0) * (N - 1.0) + 4.0 * k)) / 2.0;
    LHS = X - Si;
    RHS = GBlackScholes("p", Si, X, T, r, b, v) - (1.0 - exp((b - r) * T) * CND(-d1)) * Si / Q1;
    bi = -exp((b - r) * T) * CND(-d1) * (1.0 - 1.0 / Q1) - (1.0 + exp((b - r) * T) * ND(-d1) / (v * sqrt(T))) / Q1;
    E = 0.000001;
    // Newton Raphson algorithm for finding critical price Si
    while (fabs(LHS - RHS) / X > E)
    {
        Si = (X - RHS + bi * Si) / (1.0 + bi);
        d1 = (log(Si / X) + (b + v * v * 0.5) * T) / (v * sqrt(T));
        LHS = X - Si;
        RHS = GBlackScholes("p", Si, X, T, r, b, v) - (1.0 - exp((b - r) * T) * CND(-d1)) * Si / Q1;
        bi = -exp((b - r) * T) * CND(-d1) * (1.0 - 1.0 / Q1) - (1.0 + exp((b - r) * T) * CND(-d1) / (v * sqrt(T))) / Q1;
    }
    return Si;
}

double OptionUtil::phi(double S, double T, double gamma, double h, double i,double r, double b, double v)
{
      double lambda, kappa;
      double d;

      lambda = (-r + gamma * b + 0.5 * gamma * (gamma - 1) * v * v) * T;
      d = -(log(S / h) + (b + (gamma - 0.5) * v * v) * T) / (v * sqrt(T));
      kappa = 2.0 * b / (v * v) + 2.0 * gamma - 1.0;
      return exp(lambda) * pow(S, gamma) * (CND(d) - pow((i / S), kappa) * CND(d - 2.0 * log(i / S) / (v * sqrt(T))));

}

double OptionUtil::phi2(double S, double T2, double gamma, double h, double i, double r, double b, double v)
{
    double lambda, kappa;
    double d, d2;

    lambda = -r + gamma * b + 0.5 * gamma * (gamma - 1.0) * v * v;
    kappa = 2.0 * b / (v * v) + 2.0 * gamma - 1.0;

    d = (log(S / h) + (b + (gamma - 0.5) * v * v) * T2) / (v * sqrt(T2));
    d2 = (log(i * i / (S * h)) + (b + (gamma - 0.5) * v * v) * T2) / (v * sqrt(T2));

    return exp(lambda * T2) * pow(S, gamma) * (CND(-d) - pow((i / S), kappa) * CND(-d2));

}

double OptionUtil::ksi(double S, double T2, double gamma, double h, double I2, double I1, double t1, double r, double b, double v)
{
        double e1, e2, e3, e4;
        double f1, f2, f3, f4;
        double rho, kappa, lambda;

        e1 = (log(S / I1) + (b + (gamma - 0.5) * v * v) * t1) / (v * sqrt(t1));
        e2 = (log(I2 * I2 / (S * I1)) + (b + (gamma - 0.5) * v * v) * t1) / (v * sqrt(t1));
        e3 = (log(S / I1) - (b + (gamma - 0.5) * v * v) * t1) / (v * sqrt(t1));
        e4 = (log(I2 * I2 / (S * I1)) - (b + (gamma - 0.5) * v * v) * t1) / (v * sqrt(t1));

        f1 = (log(S / h) + (b + (gamma - 0.5) * v * v) * T2) / (v * sqrt(T2));
        f2 = (log(I2 * I2 / (S * h)) + (b + (gamma - 0.5) * v * v) * T2) / (v * sqrt(T2));
        f3 = (log(I1 * I1 / (S * h)) + (b + (gamma - 0.5) * v * v) * T2) / (v * sqrt(T2));
        f4 = (log(S * I1 * I1 / (h * I2 * I2)) + (b + (gamma - 0.5) * v * v) * T2) / (v * sqrt(T2));

        rho = sqrt(t1 / T2);
        lambda = -r + gamma * b + 0.5 * gamma * (gamma - 1.0) * v * v;
        kappa = 2.0 * b / (v * v) + (2.0 * gamma - 1.0);

        return exp(lambda * T2) * pow(S, gamma) * (CBND(-e1, -f1, rho) - pow((I2 / S), kappa) * CBND(-e2, -f2, rho)
                - pow((I1 / S), kappa) * CBND(-e3, -f3, -rho) + pow((I1 / I2), kappa) * CBND(-e4, -f4, -rho));

}

double OptionUtil::BSAmericanCallApprox2002(double S, double X, double T, double r, double b, double v)
{
    double BInfinity, B0;
    double ht1, ht2, I1, I2;
    double alfa1, alfa2, Beta, t1;

    t1 = 1.0 / 2.0 * (sqrt(5.0) - 1.0) * T;

    if (b >= r) // Never optimal to exersice before maturity
    {
        return GBlackScholes("c", S, X, T, r, b, v);
    }
    else
    {
        Beta = (1.0 / 2.0 - b / (v * v)) + sqrt(pow((b / (v * v) - 1.0 / 2.0), 2) + 2 * r / (v * v));
        BInfinity = Beta / (Beta - 1.0) * X;
        B0 = Max(X, r / (r - b) * X);

        ht1 = -(b * t1 + 2.0 * v * sqrt(t1)) * X * X / ((BInfinity - B0) * B0);
        ht2 = -(b * T + 2.0 * v * sqrt(T)) * X * X / ((BInfinity - B0) * B0);
        I1 = B0 + (BInfinity - B0) * (1.0 - exp(ht1));
        I2 = B0 + (BInfinity - B0) * (1.0 - exp(ht2));
        alfa1 = (I1 - X) * pow(I1, (-Beta));
        alfa2 = (I2 - X) * pow(I2, (-Beta));

        if (S >= I2)
        {
            return S - X;
        }
        else
        {
            return alfa2 * pow(S, Beta) - alfa2 * phi(S, t1, Beta, I2, I2, r, b, v)
                + phi(S, t1, 1, I2, I2, r, b, v) - phi(S, t1, 1, I1, I2, r, b, v)
                - X * phi(S, t1, 0, I2, I2, r, b, v) + X * phi(S, t1, 0, I1, I2, r, b, v)
                + alfa1 * phi(S, t1, Beta, I1, I2, r, b, v) - alfa1 * ksi(S, T, Beta, I1, I2, I1, t1, r, b, v)
                + ksi(S, T, 1, I1, I2, I1, t1, r, b, v) - ksi(S, T, 1, X, I2, I1, t1, r, b, v)
                - X * ksi(S, T, 0, I1, I2, I1, t1, r, b, v) + X * ksi(S, T, 0, X, I2, I1, t1, r, b, v);

        }
    }
}

double OptionUtil::BSAmericanApprox2002(string CallPutFlag, double S, double X, double T, double r, double b, double v)
{
      if(CallPutFlag=="c")
      {
              return BSAmericanCallApprox2002(S, X, T, r, b, v);
      }
      else if(CallPutFlag=="p") // Use the Bjerksund and Stensland put-call transformation
      {
              return BSAmericanCallApprox2002(X, S, T, r - b, -b, v);
      }
      else
      {
              return 0.0;
      }

}

double OptionUtil::EBSAmericanApprox2002(string OutPutFlag, string CallPutFlag, double S, double X, double T,
            double r, double b, double v, double dS)
{
      if(CallPutFlag=="p")// Then ' Value
      {
                return BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v);
      }
      else if(CallPutFlag== "d")// Then 'Delta
      {
                return (BSAmericanApprox2002(CallPutFlag, S + dS, X, T, r, b, v) - BSAmericanApprox2002(CallPutFlag, S - dS, X, T, r, b, v)) / (2.0 * dS);
      }
      else if(CallPutFlag== "e")// Then 'Elasticity
      {
                return (BSAmericanApprox2002(CallPutFlag, S + dS, X, T, r, b, v) - BSAmericanApprox2002(CallPutFlag, S - dS, X, T, r, b, v)) / (2.0 * dS) * S / BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v);
      }
      else if(CallPutFlag== "g")// Then 'Gamma
      {
                return (BSAmericanApprox2002(CallPutFlag, S + dS, X, T, r, b, v) - 2.0 * BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox2002(CallPutFlag, S - dS, X, T, r, b, v)) / (dS * dS);
      }
      else if(CallPutFlag== "gv")//Then 'DGammaDVol
      {
                return (BSAmericanApprox2002(CallPutFlag, S + dS, X, T, r, b, v + 0.01) - 2.0 * BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v + 0.01) + BSAmericanApprox2002(CallPutFlag, S - dS, X, T, r, b, v + 0.01)
                  - BSAmericanApprox2002(CallPutFlag, S + dS, X, T, r, b, v - 0.01) + 2.0 * BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v - 0.01) - BSAmericanApprox2002(CallPutFlag, S - dS, X, T, r, b, v - 0.01)) / (2.0 * 0.01 * dS * dS) / 100.0;
      }
      else if(CallPutFlag== "gp")// Then 'GammaP
      {
                return S / 100.0 * (BSAmericanApprox2002(CallPutFlag, S + dS, X, T, r, b, v) - 2.0 * BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox2002(CallPutFlag, S - dS, X, T, r, b, v)) / (dS * dS);
      }
      else if(CallPutFlag== "tg")// Then 'time Gamma
      {
                return (BSAmericanApprox2002(CallPutFlag, S, X, T + 1.0 / 256.0, r, b, v) - 2.0 * BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox2002(CallPutFlag, S, X, T - 1.0 / 256.0, r, b, v)) / (pow((1.0 / 256.0) , 2));
      }
      else if(CallPutFlag== "dddv")//Then 'DDeltaDvol
      {
                return 1.0 / (4.0 * dS * 0.01) * (BSAmericanApprox2002(CallPutFlag, S + dS, X, T, r, b, v + 0.01) - BSAmericanApprox2002(CallPutFlag, S + dS, X, T, r, b, v - 0.01)
                    - BSAmericanApprox2002(CallPutFlag, S - dS, X, T, r, b, v + 0.01) + BSAmericanApprox2002(CallPutFlag, S - dS, X, T, r, b, v - 0.01)) / 100.0;
      }
      else if(CallPutFlag== "v")// Then 'Vega
      {
                return (BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v + 0.01) - BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v - 0.01)) / 2.0;
      }
      else if(CallPutFlag== "vv")// Then 'DvegaDvol/vomma
      {
                return (BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v + 0.01) - 2.0 * BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v - 0.01)) / (0.01 * 0.01) / 10000.0;
      }
      else if(CallPutFlag== "vp")// Then 'VegaP
      {
                return v / 0.1 * (BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v + 0.01) - BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v - 0.01)) / 2.0;
      }
      else if(CallPutFlag== "dvdv")// Then 'DvegaDvol
      {
                return (BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v + 0.01) - 2.0 * BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v - 0.01));
      }
      else if(CallPutFlag== "t")// Then 'Theta
      {
                if (T <= 1.0 / 256.0)
                {
                    return BSAmericanApprox2002(CallPutFlag, S, X, 0.00001, r, b, v) - BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v);
                }
                else
                {
                    return BSAmericanApprox2002(CallPutFlag, S, X, T - 1.0 / 256.0, r, b, v) - BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v);
                }
      }
      else if(CallPutFlag== "r")// Then 'Rho
      {
                return (BSAmericanApprox2002(CallPutFlag, S, X, T, r + 0.01, b + 0.01, v) - BSAmericanApprox2002(CallPutFlag, S, X, T, r - 0.01, b - 0.01, v)) / (2.0);
      }
      else if(CallPutFlag== "fr")// Then 'Futures options rho
      {
                return (BSAmericanApprox2002(CallPutFlag, S, X, T, r + 0.01, b, v) - BSAmericanApprox2002(CallPutFlag, S, X, T, r - 0.01, b, v)) / (2.0);
      }
      else if(CallPutFlag== "f")// Then 'Rho2
      {
                return (BSAmericanApprox2002(CallPutFlag, S, X, T, r, b - 0.01, v) - BSAmericanApprox2002(CallPutFlag, S, X, T, r, b + 0.01, v)) / (2.0);
      }
      else if(CallPutFlag== "b")// Then 'Carry
      {
                return (BSAmericanApprox2002(CallPutFlag, S, X, T, r, b + 0.01, v) - BSAmericanApprox2002(CallPutFlag, S, X, T, r, b - 0.01, v)) / (2.0);
      }
      else if(CallPutFlag== "s")// Then 'Speed
      {
                return 1.0 / (dS * dS * dS) * (BSAmericanApprox2002(CallPutFlag, S + 2.0 * dS, X, T, r, b, v) - 3.0 * BSAmericanApprox2002(CallPutFlag, S + dS, X, T, r, b, v)
                              + 3.0 * BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v) - BSAmericanApprox2002(CallPutFlag, S - dS, X, T, r, b, v));
      }
      else if(CallPutFlag== "dx")// Then 'Strike Delta
      {
                return (BSAmericanApprox2002(CallPutFlag, S, X + dS, T, r, b, v) - BSAmericanApprox2002(CallPutFlag, S, X - dS, T, r, b, v)) / (2.0 * dS);
      }
      else if(CallPutFlag== "dxdx")// Then 'Strike Gamma
      {
                return (BSAmericanApprox2002(CallPutFlag, S, X + dS, T, r, b, v) - 2.0 * BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox2002(CallPutFlag, S, X - dS, T, r, b, v)) / (dS * dS);
      }
      else if(CallPutFlag== "di")//Then 'Difference in value between BS Approx and Black-Scholes Merton value
      {
                return BSAmericanApprox2002(CallPutFlag, S, X, T, r, b, v) - GBlackScholes(CallPutFlag, S, X, T, r, b, v);
      }
      else
      {
                return 0.0;
      }

}

double OptionUtil::BSAmericanCallApprox(double S, double X, double T, double r, double b, double v)
{
    double BInfinity, B0;
    double ht, i;
    double Alpha, Beta;

    if (b >= r)// Never optimal to exersice before maturity
    {
        return GBlackScholes("c", S, X, T, r, b, v);
    }
    else
    {
        Beta = (1.0 / 2.0 - b / (v * v)) + sqrt(pow((b / (v * v) - 1.0 / 2.0), 2) + 2 * r / (v * v));
        BInfinity = Beta / (Beta - 1.0) * X;
        B0 =  Max(X, r / (r - b) * X);
        ht = -(b * T + 2.0 * v * sqrt(T)) * B0 / (BInfinity - B0);
        i = B0 + (BInfinity - B0) * (1.0 - exp(ht));
        Alpha = (i - X) * pow(i, (-Beta));
        if (S >= i)
        {
            return S - X;
        }
        else
        {
            return Alpha * pow(S, Beta) - Alpha * phi(S, T, Beta, i, i, r, b, v) + phi(S, T, 1, i, i, r, b, v) - phi(S, T, 1, X, i, r, b, v) - X * phi(S, T, 0, i, i, r, b, v) + X * phi(S, T, 0, X, i, r, b, v);
        }
    }

}

double OptionUtil::BSAmericanApprox(string CallPutFlag, double S, double X, double T, double r, double b, double v)
{
      if(CallPutFlag=="c")
      {
              return BSAmericanCallApprox(S, X, T, r, b, v);
      }
      else if(CallPutFlag=="p")// Use the Bjerksund and Stensland put-call transformation
      {
              return BSAmericanCallApprox(X, S, T, r - b, -b, v);
      }
      else
      {
              return 0.0;
      }

}

double OptionUtil::EBSAmericanApprox(string OutPutFlag, string CallPutFlag, double S, double X, double T,
            double r, double b, double v, double dS)
{
      if(OutPutFlag==  "p")// Then ' Value
      {
                return BSAmericanApprox(CallPutFlag, S, X, T, r, b, v);
      }
      else if(OutPutFlag== "d")// Then 'Delta
      {
                 return (BSAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v) - BSAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v)) / (2.0 * dS);
      }
      else if(OutPutFlag== "e")// Then 'Elasticity
      {
                return (BSAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v) - BSAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v)) / (2.0 * dS) * S / BSAmericanApprox(CallPutFlag, S, X, T, r, b, v);
      }
      else if(OutPutFlag==  "g")// Then 'Gamma
      {
                return  (BSAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v) - 2.0 * BSAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag==  "gv")// Then 'DGammaDVol
      {
                return  (BSAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v + 0.01) - 2.0 * BSAmericanApprox(CallPutFlag, S, X, T, r, b, v + 0.01) + BSAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v + 0.01)
                   - BSAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v - 0.01) + 2.0 * BSAmericanApprox(CallPutFlag, S, X, T, r, b, v - 0.01) - BSAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v - 0.01)) / (2.0 * 0.01 * dS * dS) / 100.0;
      }
      else if(OutPutFlag== "gp")// Then 'GammaP
      {
                return  S / 100.0 * (BSAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v) - 2.0 * BSAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag==  "tg")// Then 'time Gamma
      {
                return  (BSAmericanApprox(CallPutFlag, S, X, T + 1.0 / 256.0, r, b, v) - 2.0 * BSAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox(CallPutFlag, S, X, T - 1.0 / 256.0, r, b, v)) / (pow((1.0 / 256.0) , 2));
      }
      else if(OutPutFlag==  "dddv")// Then 'DDeltaDvol
      {
                return  1.0 / (4.0 * dS * 0.01) * (BSAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v + 0.01) - BSAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v - 0.01)
                  - BSAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v + 0.01) + BSAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v - 0.01)) / 100.0;
      }
      else if(OutPutFlag==  "v")// Then 'Vega
      {
                return  (BSAmericanApprox(CallPutFlag, S, X, T, r, b, v + 0.01) - BSAmericanApprox(CallPutFlag, S, X, T, r, b, v - 0.01)) / 2.0;
      }
      else if(OutPutFlag==  "vv")// Then 'DvegaDvol/vomma
      {
                return  (BSAmericanApprox(CallPutFlag, S, X, T, r, b, v + 0.01) - 2.0 * BSAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox(CallPutFlag, S, X, T, r, b, v - 0.01)) / (0.01 * 0.01) / 10000.0;
      }
      else if(OutPutFlag==  "vp")// Then 'VegaP
      {
                return  v / 0.1 * (BSAmericanApprox(CallPutFlag, S, X, T, r, b, v + 0.01) - BSAmericanApprox(CallPutFlag, S, X, T, r, b, v - 0.01)) / 2.0;
      }
      else if(OutPutFlag==  "dvdv")// Then 'DvegaDvol
      {
                return (BSAmericanApprox(CallPutFlag, S, X, T, r, b, v + 0.01) - 2.0 * BSAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox(CallPutFlag, S, X, T, r, b, v - 0.01));
      }
      else if(OutPutFlag== "t" )//Then 'Theta
      {
                if (T <= 1.0 / 256.0)
                {
                        return BSAmericanApprox(CallPutFlag, S, X, 0.00001, r, b, v) - BSAmericanApprox(CallPutFlag, S, X, T, r, b, v);
                }
                else
                {
                        return BSAmericanApprox(CallPutFlag, S, X, T - 1.0 / 256.0, r, b, v) - BSAmericanApprox(CallPutFlag, S, X, T, r, b, v);
                }
      }
      else if(OutPutFlag==  "r")// Then 'Rho
      {
                return  (BSAmericanApprox(CallPutFlag, S, X, T, r + 0.01, b + 0.01, v) - BSAmericanApprox(CallPutFlag, S, X, T, r - 0.01, b - 0.01, v)) / (2.0);
      }
      else if(OutPutFlag==  "fr")// Then 'Futures options rho
      {
                return (BSAmericanApprox(CallPutFlag, S, X, T, r + 0.01, b, v) - BSAmericanApprox(CallPutFlag, S, X, T, r - 0.01, b, v)) / (2.0);
      }
      else if(OutPutFlag==  "f")// Then 'Rho2
      {
                return (BSAmericanApprox(CallPutFlag, S, X, T, r, b - 0.01, v) - BSAmericanApprox(CallPutFlag, S, X, T, r, b + 0.01, v)) / (2.0);
      }
      else if(OutPutFlag== "b")// Then 'Carry
      {
                return  (BSAmericanApprox(CallPutFlag, S, X, T, r, b + 0.01, v) - BSAmericanApprox(CallPutFlag, S, X, T, r, b - 0.01, v)) / (2.0);
      }
      else if(OutPutFlag==  "s")// Then 'Speed
      {
                return  1.0 / (dS * dS * dS) * (BSAmericanApprox(CallPutFlag, S + 2.0 * dS, X, T, r, b, v) - 3.0 * BSAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v)
                               + 3.0 * BSAmericanApprox(CallPutFlag, S, X, T, r, b, v) - BSAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v));
      }
      else if(OutPutFlag== "dx")// Then 'Strike Delta
      {
                return (BSAmericanApprox(CallPutFlag, S, X + dS, T, r, b, v) - BSAmericanApprox(CallPutFlag, S, X - dS, T, r, b, v)) / (2.0 * dS);
      }
      else if(OutPutFlag==  "dxdx")// Then 'Strike Gamma
      {
                return  (BSAmericanApprox(CallPutFlag, S, X + dS, T, r, b, v) - 2.0 * BSAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BSAmericanApprox(CallPutFlag, S, X - dS, T, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag==  "di")// Then 'Difference in value between BS Approx and Black-Scholes Merton value
      {
                return BSAmericanApprox(CallPutFlag, S, X, T, r, b, v) - GBlackScholes(CallPutFlag, S, X, T, r, b, v);
      }
      else if(OutPutFlag==  "BSIVol")// Then 'Equivalent Black-Scholes-Merton implied volatility
      {
                if (S >= S * exp(b * T))
                {
                    return GBlackScholesImpVolBisection("c",S,X,T,r,b,BSAmericanApprox(CallPutFlag, S, X, T, r, b, v));
                }
                else
                {
                    return GBlackScholesImpVolBisection("p",S,X,T,r,b,BSAmericanApprox(CallPutFlag, S, X, T, r, b, v));

                }
      }
      else
      {
                return 0.0;
      }

}

double OptionUtil::BAWAmericanCallApprox(double S, double X, double T, double r, double b, double v)
{
    double Sk, N, k;
    double d1, Q2, a2;

    if (b >= r)
    {
        return GBlackScholes("c", S, X, T, r, b, v);
    }
    else
    {
        Sk = Kc(X, T, r, b, v);
        N = 2.0 * b / (v * v);
        k = 2.0 * r / (v * v * (1 - exp(-r * T)));
        d1 = (log(Sk / X) + (b + v * v * 0.5) * T) / (v * sqrt(T));
        Q2 = (-(N - 1.0) + sqrt((N - 1.0) * (N - 1.0) + 4.0 * k)) / 2.0;
        a2 = (Sk / Q2) * (1.0 - exp((b - r) * T) * CND(d1));
        if (S < Sk)
        {
            return GBlackScholes("c", S, X, T, r, b, v) + a2 * pow((S / Sk), Q2);
        }
        else
        {
            return S - X;
        }
    }

}

double OptionUtil::BAWAmericanPutApprox(double S, double X, double T, double r, double b, double v)
{
    double Sk, N, k;
    double d1, Q1, a1;

    Sk = Kp(X, T, r, b, v);
    N = 2 * b / (v * v);
    k = 2 * r / (v * v * (1 - exp(-r * T)));
    d1 = (log(Sk / X) + (b + v * v * 0.5) * T) / (v * sqrt(T));
    Q1 = (-(N - 1.0) - sqrt((N - 1.0) * (N - 1.0) + 4.0 * k)) / 2.0;
    a1 = -(Sk / Q1) * (1.0 - exp((b - r) * T) * CND(-d1));

    if (S > Sk)
    {
        return GBlackScholes("p", S, X, T, r, b, v) + a1 * pow((S / Sk), Q1);
    }
    else
    {
        return X - S;
    }

}

double OptionUtil::BAWAmericanApprox(string CallPutFlag, double S, double X, double T, double r, double b, double v)
{
      if(CallPutFlag=="c")
      {
              return BAWAmericanCallApprox(S, X, T, r, b, v);
      }
      else if(CallPutFlag=="p")
      {
              return BAWAmericanPutApprox(S, X, T, r, b, v);
      }
      else
      {
              return 0.0;
      }

}

double OptionUtil::EBAWAmericanApprox(string OutPutFlag, string CallPutFlag, double S, double X, double T,
            double r, double b, double v, double dS)
{
      if(OutPutFlag== "p")// Then ' Value
      {
                return BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v);
      }
      else if(OutPutFlag== "d")// Then 'Delta
      {
                return (BAWAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v) - BAWAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v)) / (2.0 * dS);
      }
      else if(OutPutFlag== "e")// Then 'Elasticity
      {
                return (BAWAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v) - BAWAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v)) / (2.0 * dS) * S / BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v);
      }
      else if(OutPutFlag== "g")// Then 'Gamma
      {
                return (BAWAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v) - 2.0 * BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BAWAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag== "gv")// Then 'DGammaDVol
      {
                return (BAWAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v + 0.01) - 2.0 * BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v + 0.01) + BAWAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v + 0.01)
                          - BAWAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v - 0.01) + 2.0 * BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v - 0.01) - BAWAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v - 0.01)) / (2.0 * 0.01 * dS * dS) / 100.0;
      }
      else if(OutPutFlag== "gp")// Then 'GammaP
      {
                return S / 100.0 * (BAWAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v) - 2.0 * BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BAWAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag== "tg")// Then 'time Gamma
      {
                return (BAWAmericanApprox(CallPutFlag, S, X, T + 1.0 / 256.0, r, b, v) - 2.0 * BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BAWAmericanApprox(CallPutFlag, S, X, T - 1.0 / 256.0, r, b, v)) / (pow((1.0 / 256.0), 2));
      }
      else if(OutPutFlag== "dddv")// Then 'DDeltaDvol
      {
                return 1.0 / (4.0 * dS * 0.01) * (BAWAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v + 0.01) - BAWAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v - 0.01)
                           - BAWAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v + 0.01) + BAWAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v - 0.01)) / 100.0;
      }
      else if(OutPutFlag== "v")// Then 'Vega
      {
                return (BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v + 0.01) - BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v - 0.01)) / 2.0;
      }
      else if(OutPutFlag== "vv")// Then 'DvegaDvol/vomma
      {
                return (BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v + 0.01) - 2.0 * BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v - 0.01)) / (0.01 * 0.01) / 10000.0;
      }
      else if(OutPutFlag== "vp")// Then 'VegaP
      {
                return v / 0.1 * (BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v + 0.01) - BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v - 0.01)) / 2.0;
      }
      else if(OutPutFlag== "dvdv")// Then 'DvegaDvol
      {
                return (BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v + 0.01) - 2.0 * BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v - 0.01));
      }
      else if(OutPutFlag== "t")// Then 'Theta
      {
                if (T <= 1.0 / 256.0)
                {
                    return BAWAmericanApprox(CallPutFlag, S, X, 0.00001, r, b, v) - BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v);
                }
                else
                {
                    return BAWAmericanApprox(CallPutFlag, S, X, T - 1.0 / 256.0, r, b, v) - BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v);
                }
      }
      else if(OutPutFlag== "r")// Then 'Rho
      {
                return (BAWAmericanApprox(CallPutFlag, S, X, T, r + 0.01, b + 0.01, v) - BAWAmericanApprox(CallPutFlag, S, X, T, r - 0.01, b - 0.01, v)) / (2.0);
      }
      else if(OutPutFlag== "fr")// Then 'Futures options rho
      {
                return (BAWAmericanApprox(CallPutFlag, S, X, T, r + 0.01, b, v) - BAWAmericanApprox(CallPutFlag, S, X, T, r - 0.01, b, v)) / (2.0);
      }
      else if(OutPutFlag== "f")// Then 'Rho2
      {
                return (BAWAmericanApprox(CallPutFlag, S, X, T, r, b - 0.01, v) - BAWAmericanApprox(CallPutFlag, S, X, T, r, b + 0.01, v)) / (2.0);
      }
      else if(OutPutFlag== "b")// Then 'Carry
      {
                return (BAWAmericanApprox(CallPutFlag, S, X, T, r, b + 0.01, v) - BAWAmericanApprox(CallPutFlag, S, X, T, r, b - 0.01, v)) / (2.0);
      }
      else if(OutPutFlag== "s")// Then 'Speed
      {
                return 1.0 / (dS * dS * dS) * (BAWAmericanApprox(CallPutFlag, S + 2.0 * dS, X, T, r, b, v) - 3.0 * BAWAmericanApprox(CallPutFlag, S + dS, X, T, r, b, v)
                             + 3.0 * BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v) - BAWAmericanApprox(CallPutFlag, S - dS, X, T, r, b, v));
      }
      else if(OutPutFlag== "dx")//Then 'Strike Delta
      {
                return (BAWAmericanApprox(CallPutFlag, S, X + dS, T, r, b, v) - BAWAmericanApprox(CallPutFlag, S, X - dS, T, r, b, v)) / (2.0 * dS);
      }
      else if(OutPutFlag== "dxdx")//Then 'Strike Gamma
      {
                return (BAWAmericanApprox(CallPutFlag, S, X + dS, T, r, b, v) - 2.0 * BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v) + BAWAmericanApprox(CallPutFlag, S, X - dS, T, r, b, v)) / (dS * dS);
      }
      else if(OutPutFlag== "di")// Then 'Difference in value between BS Approx and Black-Scholes Merton value
      {
                return BAWAmericanApprox(CallPutFlag, S, X, T, r, b, v) - GBlackScholes(CallPutFlag, S, X, T, r, b, v);
      }
      else
      {
                return 0.0;

      }
}

double OptionUtil::BAWAmericanImpVolBisection(string CallPutFlag, double S, double x, double T, double r, double b, double cm)
{
    double vLow, vHigh, vi;
    double cLow, cHigh, epsilon;
    int counter;

    vLow = 0.005;
    vHigh = 4;
  //  epsilon = 0.00000001;
    epsilon = 0.0001;

    cLow = BAWAmericanApprox(CallPutFlag, S, x, T, r, b, vLow);
    cHigh = BAWAmericanApprox(CallPutFlag, S, x, T, r, b, vHigh);
    counter = 0;
    vi = vLow + (cm - cLow) * (vHigh - vLow) / (cHigh - cLow);
    while (fabs(cm - BAWAmericanApprox(CallPutFlag, S, x, T, r, b, vi)) > epsilon)
    {
        counter = counter + 1;
        if (counter == 500) return 0.0;
        if (BAWAmericanApprox(CallPutFlag, S, x, T, r, b, vi) < cm)
        {
            vLow = vi;
        }
        else
        {
            vHigh = vi;
        }
        cLow = BAWAmericanApprox(CallPutFlag, S, x, T, r, b, vLow);
        cHigh = BAWAmericanApprox(CallPutFlag, S, x, T, r, b, vHigh);
        vi = vLow + (cm - cLow) * (vHigh - vLow) / (cHigh - cLow);
    }
    return vi;
}

double OptionUtil::BSAmerican2002ImpVolBisection(string CallPutFlag, double S, double x, double T, double r, double b, double cm)
{
    double vLow, vHigh, vi;
    double cLow, cHigh, epsilon;
    int counter;

    vLow = 0.005;
    vHigh = 4;
    //  epsilon = 0.00000001;
    epsilon = 0.0001;

    cLow = BSAmericanApprox2002(CallPutFlag, S, x, T, r, b, vLow);
    cHigh = BSAmericanApprox2002(CallPutFlag, S, x, T, r, b, vHigh);
    counter = 0;
    vi = vLow + (cm - cLow) * (vHigh - vLow) / (cHigh - cLow);
    while (fabs(cm - BSAmericanApprox2002(CallPutFlag, S, x, T, r, b, vi)) > epsilon)
    {
        counter = counter + 1;
        if (counter == 500)
        {
            if (fabs(cm - BSAmericanApprox2002(CallPutFlag, S, x, T, r, b, vi)) < 0.01)
            {
                return vi;
            }
            else
            {
                return 0.0;
            }
        }
        if (BSAmericanApprox2002(CallPutFlag, S, x, T, r, b, vi) < cm)
        {
            vLow = vi;
        }
        else
        {
            vHigh = vi;
        }
        cLow = BSAmericanApprox2002(CallPutFlag, S, x, T, r, b, vLow);
        cHigh = BSAmericanApprox2002(CallPutFlag, S, x, T, r, b, vHigh);
        vi = vLow + (cm - cLow) * (vHigh - vLow) / (cHigh - cLow);
    }
    return vi;
}